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ABSTRACT 



The current antenna used in the Marine Radio AN/MRC-142 is a focal point fed 
reflector with a first sidelobe of -10 dB. A new antenna design is required with a sidelobe 
reduction of another 10 dB, and a minimum impact on the system upgrade cost. An 
improved feed consisting of a fed dipole and disk in the vicinity of a parasitic ring was 
designed for use with the current reflector. The disk acts as a ground plane for the dipole, 
with the ring narrowing the beamwidth, resulting in a stronger edge taper on the reflector. 

An existing computer code based on the the method of moments was used to 
compute the gain and the radiation pattern. Circular parabolic reflectors of 36 inches and 
48 inches were analyzed. The results show that the 36 inch reflector cannot satisfy the 
conditions over the entire frequency band (1350 to 1850 MHz) simultaneously. The 48 
inch reflector came very close to satisfying all the conditions, although the sidelobe level 
is off by approximately 2 dB at wide angles in the H-plane. 
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INTRODUCTION 



A. BACKGROUND 

An antenna is the transducer that intercepts propagating electromagnetic waves and 
directs them into guided waves in a transmission line. In the case of a transmitting 
antenna, guided electromagnetic energy in a transmission line is converted into radiated 
electromagnetic waves. In short, the antenna couples the system to the environment [Ref. 
!]• 

Almost all radar, communication, and electronic warfare (EW) systems require 
some form of antenna. Even for highly directive antennas, all of the energy is not 
concentrated into the main beam; there is always radiation in directions other than the 
intended one. These sidelobes, as they are called, are of particular importance in 
electronic warfare, since they represent a major limitation to the system's effectiveness. 
Jamming signals or false targets received through a sidelobe can cause confusion and 
failure of a mission more easily than false signals detected by the main beam. The 
reduction and control of the position and magnitude of sidelobes is extremely important 
when attempting to localize target azimuths in a hostile electronic environment. 

Antenna sidelobe control is one of the most basic methods of defeating jammers 
and avoiding detection. Low sidelobes also reduce multipath signals, which can cause 
signal fading in communications systems and tracking errors in radars. The methods of 
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sidelobe suppression for passive antennas are well known, and the technology is low risk. 
However, low sidelobes are accompanied by a reduction in gain, and generally require a 
more complicated antenna design. Consequently, there are tradeoffs in the antenna design 
that impact the overall system performance and cost. 

B. OBJECTIVE 

This thesis investigates the feasibility of reducing antenna sidelobes for the Marine 
Radio MRC-142. This system provides medium-range ultra high frequency (UHF) 
line-of-sight multi-channel radio communications over a distance of up to 35 miles. The 
AN/MRC-142 radio terminal set is installed in the M-998 variant of the High Mobility 
Multipurpose Wheeled Vehicle (HMMWV) and when deployed appears as in Figure 1. 

This unit can be configured as a repeater to extend the communications range 
between two other systems or to overcome terrain conditions that prevent a direct 
line-of-sight path between the other two pieces of equipment. The AN/MRC-142 
interfaces with a distant end AN/MRC-142 over a 1350 to 1850 MHz UHF radio link 
using parabolic antennas that are mounted on masts which can be extended up to 50 feet. 
The current design uses a parabolic reflector with a focal-point feed and has a first 
sidelobe requirement of -10 dB relative to the mainbeam. A reduction of another 10 dB is 
desired for the new design, with a minimum impact on the system upgrade cost. 
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50 Foot Telescoping Antenna 




Figure 1. Deployed AN/MRC-142 Radio Terminal Set [Ref. 2]. 



C. METHOD OF SOLUTION 

Possible options for the low sidelobe antenna include: 

1 . an array antenna; 

2. a completely redesigned reflector and feed, or; 

3. an improved feed used in conjunction with the current reflector. 

The complete antenna specifications are given in the Appendix. Options 1 and 2 could 
satisfy the specifications in their entirety, but the cost of design and development would 
be significant. The third option is the one pursued here. Even though all of the antenna 
specifications are not achieved simultaneously with the new feed, there is significant 
improvement in all areas. Most importantly, the cost of the new feed is minimal because 
of its simplicity. 
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The proposed feed design consists of a fed dipole in the vicinity of a parasitic ring 
and disk as shown in Figure 2 [Ref. 3], Essentially the disk acts as a ground plane for the 
dipole and the ring narrows the beamwidth, resulting in a stronger edge taper on the 
reflector. The dimensions are chosen to minimize the sidelobes. This "tuning" approach 
is narrow band and, as discussed in Chapter III, does not allow an optimum feed design 
for both the transmit and receive frequencies. 




Figure 2. Dipole-Disk Antenna with Ring Used as a Feed in Paraboidal Reflector. 



D. SCOPE OF THESIS 



In Chapter II the method of analysis is discussed in detail, 
pattern data presented are generated using the method of moments. 



All of the gain and 
Chapter III discusses 



4 



the tradeoffs involved in the antenna design. As mentioned earlier, the greatest challenge 
was selecting the geometrical parameters to provide acceptable performance at both 
frequencies (1350 MHz and 1850 MHz). Finally, Chapter IV presents a summary of the 
work performed along with some conclusions. 
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II. THEORETICAL BACKGROUND 



A. E-FIELD INTEGRAL EQUATION AND THE METHOD OF MOMENTS 

To determine the pattern and gain of a reflector antenna, it is first necessary to find 
the current on the reflector surfaces. This is accomplished by solving an integral or 
differential equation for the current which is derived from Maxwell's equations and the 
boundary conditions. For wire and surface geometries of arbitrary complexity, the 
integral or differential equation must be solved numerically on a high speed digital 
computer. The method of moments (MM) is such a technique which reduces the E-field 
integral equation (EFIE) to a matrix equation which is then solved on a computer [Ref. 4], 
The following paragraphs summarize the method of moments technique. Phaser notation 
is used throughout and e 1 ” 1 time convention is assumed and suppressed. 

If the current on a body is known, then the scattered (or radiated) field can be 
determined from [Ref. 5] 



where co = 27tf is the radian frequency, andp and e are the permeability and permittivity 
of free space. A is the magnetic vector potential, which can be represented as 



E = -V4> - jffl A = -jco A - j^V(V • A), 



( 1 ) 



A = 




( 2 ) 



where J s is the surface current and S refers to the surface of the conductor. 
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The electric field integral equation is based on the boundary' condition that the total 
tangential electric field on a perfectly electric conducting (PEC) surface of an antenna or 
scatterer is zero. This can be expressed as 

E}(r = r s ) = Ej(r = r s ) + E, s (r = r s ) = 0 on S (3) 

where S is the conducting surface of the antenna or scatterer, as stated before, andr = r s 
is the distance from the origin to any point on the surface of the antenna or scatterer. The 

subscript "t" indicates tangential components. 

Referring to Figure 3, Equations (1) and (2) can be combined and expressed as 





Figure 3. Coordinate Systems for Computing Radiation Fields. 
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E‘(F) = -jJ[P I IJ,J,(F')G(r,r')ds' + Vl| 1 V'.j l (f')G(F,r')ds'J, 



W 



where 



GO 7 ,?') 



e -J kR 
4nR ’ 



( 5 ) 



and 



R=|r-r'l- (6) 

Now, applying the boundary condition of Equation (3) gives 

E s (r) = E t s (r s ) = J j[3 2 ll s J s (7 / )G(7,r / )ds / + Vll s V / .J s (r')G(7,7 / )ds / ] (7) 

This is an integral equation with J s as the unknown Equation (7) is referred to as the 
electric-field integral equation. In Equation (7),V and V' are the gradients with respect to 

the observation (unprimed) and source (primed) coordinates, respectively, and 
G( 7,7 ) 

is the Green's function for a three dimensional radiator or scatterer 

For convenience, an operator L can be defined which will operate on 1, 

L(J S ) =jcopJ J s J S G; ds' + £v[vJj s T s G; ds' +JJ S I S *VG; ds' (8) 

where unknown current J s depends only on the primed quantities. G is a scalar quantity 
that depends on both the primed and unprimed quantities. It is desirable to have all the 

differential operations in terms of the primed variables of integration. Equation (8) can 

therefore be written as 

L(J S ) =jop|J s J S G, ds' + £v[jJ s G; V. T s ds' + JJ s J s • V; G; ds', (9) 
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where the vector identity 



V(J s • VG) = (J s • V)VG + (VG • V)J s + J s x (V x VG) + VG x (V x J s ) . (10) 

has been used. 

The gradient of G in terms of the two sets of coordinates is related by 

VG = -0k + })^? = -V'G (11) 

Using the surface divergence theorem developed by Mautz and Harrington [Ref. 5], it can 
be shown that 

J7 S V' • (J s G)ds / = 0 (12) 

Finally, Equation (7) can be written as 

L(J S ) = j J s Ljor|JsG - £V[(V' . J s )G]]ds'. (13) 

Now the electric field integral equation can be represented as 

Ej(r s ) = L(I s ). (14) 

Equation (14) can be solved numerically using the method of moments The current 
on the antennas is expanded in terms of a finite series of known basis functions with 
unknown amplitudes. The antenna under consideration consists of surfaces and thin 
wires, so basis functions must be defined on each. On a surface two orthogonal 
components are required. The coordinate system for the reflector is shown in Figure 4 
and the following basis functions defined on the surface [Ref. 6]: 
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SEGMENTED REFLECTOR GENERATING CURVE 




Figure 4. A Reflector Geometry for Method of Moments. 



n = 0,±1, ... ±oc, i = 1, 2, . . ,,N S - 2 
j = n = 0,±1, ... ±oo, i = 1,2, ...,N S - 1 



(15) 



T,(t) is the triangle function (which extends over two segments, i and i + 1) and P,(t) 
is the pulse function A point on the surface of the antenna is specified by the coordinates 
(t, 4>), where t is an arclength variable along the reflector generating curve. The distance 
of a point from the axis of symmetry (z-axis) is given by p and N s is the number of 
surface generating points. A discussion leading to the choice of these basis functions is 
given by Mautz and Harrington [Ref. 5], 

For thin wires the radius "a" is much less than a wavelength and the current will 

A 

only have an axial component. Thus the operator V becomes t(d/<5t), where t is the arc 
length along the wire. Furthermore, the current is constant around the wire and J/ S ds' 
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reduces to 27taJ,dt\ as shown in Figure 5. Thus appropriate basis functions on the wire 
are [Ref. 7] 




Figure 5. Thin Wire Approximation. 



Ji = t 



T,(t) 

2na ’ 



i=l,2...,N w . 



(16) 



A 

In this case t is a unit vector along the wire axis and H is the number of wire generating 
points. The total current is the sum over all expansion functions 



N w -2 - oc 

T _ y twtw . y 
J S — ^ 1 m J m 

m=l n — x 



N s -2 _ Ns-1 a - 6 

S Ip n Jpn + ^ Iqn J qn 
p=l V F q=l 



(17) 



The coefficients I are the unknowns which need to be determined 

From the linearity of L, Equation (14) can be combined with Equation (17) as 
follows 



E{= SliL(Jj). 

1=1 
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(18) 



where N b is the total number of basis functions of both types. N b equations need to be 

provided to enable a solution for the unknown constants I,. Galerkin’s method is used to 

reduce Equation (18) to a matrix equation, and weighting functions are defined such that 
Wk = J k (Ref- 8, 9], Each side is multiplied by each of the basis functions and then each 
side integrated. The result is N b equations of the form 

JJ, kl E i .W k ds = [JJ Sk > k M 1 L(j 1 )+J| si2 W k .I 2 L(7 2 ) + ...+Jl SkN W k .I N L(7 N )]ds 



k=l,2,...N b (19) 

which can be written as 

1 J E‘ • W k ds = JJ [W k • I I,L(7 ,)ds] (20) 

-K M 1=1 

Interchanging summation and integration, 

j L E' .W k ds = i b lJJW k . L(7 i)ds . (21) 

k i=l 

Equation (21) can be cast in the form 

V k =IZ ik I,. (22) 

i=l 

for k = 1,2, ...,N b , where the impedance elements are given by 

Z lk = JJ S dsJJ s [jcopWk • 7i - i(V ' • 7 i)(V . W k )]Gds ' . (23) 

The excitation elements are defined by 

V k =H s W k .E i ds, (24) 
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In matrix notation Equation (22) becomes 



M = [Z][I], (25) 

and the unknown coefficients are determined by solving the matrix equation 

[1] = [Z] '[V] (26) 

. . • A A 

Submatrices in Equation (25) associated with the t and <f> surface components can be 
identified as follows 



" [V'l ‘ 


= 


[z;hz?i " 


In " 


[vt] 


- 


[zS’](z7] 


I* 

In 



(27) 



The impedance matrix will have the following block structure when all of the 
surface and wire basis functions are considered 



7 ss 

a -m 




Z ws rj ws y \vs 

-M 



ySW 




" I-M ' 




1 

Crt 1 

> 

1 


ySW 




n 


= 


Vo 


7SU 




III 




vj, 


7 \V\Y 


L i w J 




V' v 



(28) 



The superscripts s and w refer to surface and wire, respectively. Although n in 
Equation (17) can take on values between! oc, it has been truncated at! M in the above 
matrix. In practice, when the feed is on the reflector axis of symmetry, a converged 
solution can be obtained with M =1. If the pattern is calculated for the antenna 
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transmitting, the excitation vector will have only one non-zero value corresponding to the 
dipole feed point. 

The unknown expansion coefficients are determined by solving the matrix in 
Equation (27) for [I], To compute the electric field, the coefficients are used in Equation 
(17), which in turn is used in Equations (1) and (2). Define a measurement vector, 



where E r denotes a unit radiated plane wave weighted by the current on S. In general J s 
will give rise to both 0 and 4> components in the far field. The radiated 0 component is of 

the form 



where u = sin0cos<}>, v = sin0sins(J>, and w = cos0. By substituting Equation (30) into 
Equation (29), there results 



The total 0 component of the electric field is the superposition of all surface current 
contributions 



Rm =JJ S E r • J s ds, 



(29) 



E r = 0e~ jk,r = 0e“ jk|xu+yv+2u| , 



(30) 



R® =JJ S e- jk|xu+yv+zw) (0*J s )ds . 



(31) 




(32) 



The analysis holds true for the <j> component. The measurement elements are 



Rm = j j s e“j k, ™ +yv+zw| (<t> • J S )ds , 



(33) 
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and the electric field is then defined by 



E * = : S !e ' ik, 5 R ” 1 " 



( 34 ) 



B. GAIN CALCULATION 

The degree to which any particular antenna pattern is concentrated into a beam is 
known as directivity' which is defined by 



pi P max 

i) > 

t ave 



(35) 



where P max is the maximum power flux radiated, and P ave is the average power flux 
radiated Since the average power flux radiated is equal to the total power radiated (P, ot ) 
divided by 4n steradians, 

D = 47t—v (36) 

Directivity, as defined above, is a function of the antenna pattern only, and does not take 
into account the antenna efficiency or losses. In order to compare the relative 
performance of antennas in any electronic-warfare system, it is necessary to consider 
antenna efficiency. The pertinent quantity is gain G which is defined as 



G = KD, 



(37) 
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where K is the efficiency factor (0<K< 1). Thus it is seen that for a 100% efficient 
antenna (K = 1), the gain would be equal to the directivity. For all practical antennas K < 
1 . Gain is frequently expressed as a decibel ratio 

G,db = 101og ]0 G (38) 

A very important relationship exists between gain, physical antenna dimensions, 
and electrical wavelength, expressed by 

G = 47t (39) 

where A is the "effective" area of the antenna and^. is the wavelength. The effective area 
of the antenna and the wavelength must be expressed in the same units (meters, inches, 
feet, etc ). The total radiated power required in Equation (36) can be obtained by 
integrating the electric field over a sphere 

P,ot = ? j { |Ee | 2 + |E*| 2 }r 2 sin 0d0d<|> (40) 

o o 

where E e and E^, are known from Equations (32) and (34). 
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III. LOW SIDELOBE ANTENNA DESIGN STUDY 



A. SUMMARY OF TRADEOFFS IN THE ANTENNA DESIGN 

The purpose of this thesis is to achieve a reduction in the antenna sidelobes for the 
Marine radio AN/MRC-142. The approach will be to upgrade the existing reflector in 
Figure 6 by implementing a feed design similar to that of Kildal and Skyttemyr [Ref. 3] 
shown in Figure 2 of Chapter I. 




Figure 6. AN/MRC-142 Antenna Assembly and Offset Adapter [Ref. 2.]. 

The size of the reflector aperture is 36 x 48 inches as shown in Figure 7. The 
projected aperture is not circular, but has two circular segments removed The method of 
moments code that was used to compute the gain and radiation pattern can only handle 
rotationally symmetric reflectors. In the design of the antenna, circular parabolic 
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reflectors of 36 and 48 inches were analyzed. Additionally, the structure supports only 
vertical currents, whereas the method of moments code assumes radial and azimuthal 
currents (solid surface). The actual AN/MRC-142 performance will lie somewhere 
between the two. The pattern will resemble the 48 inch reflector in the plane of the wide 
dimension, and the 36 inch reflector in the plane of the narrow dimension. The effect on 
gain can be accurately estimated based on the area reduction from the 48 inch diameter 
circular antenna due to removal of the circular segments. 



Figure 7. Antenna Radius Determination. 

The reflector surface is a section of a parabaloid and is fed from the front. Four 

factors affect the efficiency of this type of antenna: 

1. Spillover efficiency; 

2. Illumination efficiency; 

3. Blockage by the feed; 




h — rr — H 
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4. Cable and line losses from the antenna input to the feed. 



Taking into consideration all these factors, most of the burden in maximizing the 
efficiency is placed upon the antenna feed. The dipole-disk antenna, which essentially 
consists of a halfwave dipole over a circular plane reflector, is widely used as feed for 
paraboloidal reflector antennas at VHF and UHF. A method of improving this feed 
consists in incorporating a circular conducting beam-forming ring (BFR) over the dipole 
in a plane parallel to the disk The beam-forming ring compresses the H-plane pattern of 
the dipole-disk feed, whereas it has no significant effect on the E-plane pattern When 
this antenna is used as feed for paraboloidal reflectors, it exhibits near-identical principal 
plane aperture illumination and reduces H-plane spillover lobes. 

The diameter of the disk should be as small as possible in order to reduce blockage 
loss. However a small disk and ring result is a broad feed radiation pattern which does 
not provide sufficient edge taper for a low sidelobe pattern. To obtain the best 
performance a parametric study was performed for the geometrical quantities shown in 
Figure 8 and listed in Table 1. The parameters and their effect in antenna radiation 
pattern and gain discussed in the following sections. All dimensions used by the 
computer code are in wavelengths. 

1. Reflector Size 

The size of the reflector is 36 x 48 inches as shown in Figure 7. Two different 
designs were obtained for each frequency: one for 36 inches and the other one for 48 
inches. At 1350 MHz, the wavelength is 8.75 inches and at 1850 MHz the wavelegth is 
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6.38 inches. Thus 36 inches corresponds to 4.1 A. and 5.6X, while 48 inches corresponds 
to 5.5A. and 1.5X at these two frequencies. 

The focal-length-to-diameter ratio (f/D) of the existing reflector antenna could 
not be determined from AN/MRC-142 technical manuals. Typical values are in the range 
of 0.3 to 0.55. Smaller values are preferable to reduce spillover loss and increase edge 
taper. A value of 0.4 was chosen for f/D. This parameter could be varied slightly without 
significant modification of the antenna's gain and pattern. 




Figure 8. Program Parameters 



20 



DM 

DS 

ZC 

FOD 

ZR 

RC 

SW 

DLN 

PLN 

DDIP 

AW 

HLGTH 

DHP 



TABLE 1. DESIGN PARAMETERS 

diameter of paraboloid 
diameter of disk 

distance of dipole to ground plane 
focal length to diameter ratio of the parabolid 
distance of ring from fed dipole along the z axis 
radius of ring 

ring width (modeled as a thin strip) 

fed dipole length 

parasitic director dipole length 

distance between DLN and PLN along the z axis 

radius of dipoles 

parasitic (H plane) dipole length (not used) 

H plane dipole distance from feed along the y axis (not used) 



2. Gain Calculation 

Gain estimates for the reflector and feed were obtained by first integrating the 
total electric field to find the total radiated power P tot . The result is then used in Equation 

[36] to determine the directivities shown in Table 2. These numbers include spillover and 
aperture efficiency because the total field is integrated (that is, feed plus scattered fields). 

Losses that are not accounted for and thus determine the value of K in Equation 

[37] include: 

1. Transmission line losses for the cable that feeds the dipole. Assuming 0.5 dB/ft 
and 3 ft of cable yields 1.5 dB of loss. 
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2. Reflection loss from the reflector surface. This is due to the fact that the surface is 
a grid rather than a solid. A conservative estimate is 0.5 dB. 

3. Any impedence mismatch loss. A VSWR of 2:1 results is approximately 0.5 dB 
of loss. 

Based on the above numbers the total loss is about 2.5 dB. The corresponding 
gains are shown in Table 3. 



TABLE 2 


. CALCULATED DIRECTIVITIES. 




48 inches 


36 inches 


1350 MHz 


24.75 dB 


22.19 dB 


1850 MHz 


27.44 dB 


24.90 dB 




TABLE 3. 


CALCULATED GAINS WITH LOSS. 




48 inches 


36 inches 


1350 MHz 


22.25 dB 


19.19 dB 


1850 MHz 


23.94 dB 


21.90 dB 



B. MODIFICATIONS TO ANTENNA DESIGN 

Each of the parameters of the feed and reflector were varied in small increments and 
their effect on the antenna pattern and gain was evaluated The results are summarized 
below. 



22 



1. Baseline Feed Design 

The dimensions of the baseline feed design are listed in Table 4 The disk serves 
as a ground plane for the dipole. The phase center of the feed will vary with frequency 
and angle. If the disk is an effective ground plane the phase center of the dipole and its 
image will be approximately at the midpoint of the two. However in the following 
discussion the dipole is considered to be the center of the feed. 

Low sidelobes require a strong amplitude taper across the reflector, which in turn 
requires a narrow feed beamwidth. Generally a narrow feed beamwidth is associated with 
a feed aperture which results in increased blockage. Thus there is a tradeoff between 
sidelobe taper and blockage loss. 

The geometrical parameters in Table 4 were the starting point for the parametric 
study, which is summarized in the following sections The baseline antenna performance 
is shown in Figure 10. The solid line is the E-plane (4>=0°) pattern and the dashed line the 



H-plane (<j)=90 o ) pattern. 




0, degrees 



Figure 10. Baseline Antenna Pattern. 



TABLE 4. BASELINE FEED GEOMETRY 

Z distance or 

Variable Value coordinate 1 


disk diameter 


4.1 


Zc = .25 


ring diameter 


0.7 


Zr = 0.0 


fed dipole length 


0.5 


Zd = -0.25 


parasitic dipole length 


not part of baseline design 


fed dipole height above disk 


0.5 


Zc-Zd 


ring heigth above disk 


0.25 


Zc-Zr 


wire radius 


.OOU 


- 


1 The fed dipole is located at Z = 0; positive Z is away from the reflector 
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2. Disk Size and Location 



The disk must be large enough to serve as an effective ground plane for the 
dipole. Generally the dipole height must be significantly less than the diameter of the 
disk. For the dipole and its image to add constructively in the direction of the reflector, 
the spacing should be approximately .25X. The patterns for disk locations of Z = +. IX, 
and -.IX are shown in Figures 11 and 12. All other feed parameters are held constant 
(values given in Table 4). The reflector diameter is 4. IX and f/D = 0.5 in all cases. The 
high levels in the H plane patterns for 80° < 0 < 120° are due to feed spillover. A 
comparison of the figures shows that there is a tradeoff between spillover levels in the 
rear hemisphere (90 < 0 < 180) and sidelobe levels in the forward hemisphere (0< 90°). 
The impact of changing the disk radius by 1 X is illustrated in Figures 13 and 14. 
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0 




Figure 11. Radiation Pattern with Disk Shift of +0.1 X. 
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Figure 12. Radiation Pattern with Disk Shift of -0.1 X. 
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Figure 13. Radiation Pattern with Disk Size .IX.. 
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0, degrees 



Figure 14. Radiation Pattern with Disk Reduced -.IX. 
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3. Ring Size and Position 

The beam-forming ring (BFR) compresses the H-plane beamwidth but leaves the 
E-plane beamwidth essentially unaffected. Thus the ring size and location can be used to 
equalize the feed beamwidth. Patterns for two ring radii are shown in Figures 15 and 16 
The ring acts somewhat like a director in a Yagi array when it is located in front 
of the dipole (towards the reflector) [Ref. 6], When the ring is located behind the dipole it 
has little effect on the reflector scattered field. As illustrated in Figure 17 and 18, a shift 
of .R results in about a 1 dB difference in the first sidelobe level. 
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Figure 15. Radiation Pattern with Ring Size of 0.6X. 
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Figure 16. Radiation Pattern with Ring Size of 0.5\. 
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Figure 17. Radiation Pattern with Ring at Z = -0.2X. 
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0 




Figure 18. Radiation Pattern with Ring at Z = -0.1 A.. 

4. Adding a Parasitic Dipole 

A parasitic dipole was added as a director in order to improve the directivity of 
the feed and therefore increase the reflector edge taper. This method is effective in 
narrowing the beamwidth of cavity-backed dipole antennas. In this case it was found that 
changing the size and position of the parasitic element has no significant effect on the 
performance of the antenna and does not offer any advantage. This is probably because 
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the parasitic dipole is not located within a cavity. Adding a second ring around the 
parasitic dipole could improve the performance. In an effort to keep the feed as simple as 
possible, the parasitic element was not considered further. 

5. Defocusing the Feed 

The phase center of the feed varies with frequency and angle, but is generally 
located between the fed dipole and disk. The antenna bandwidth can be increased by 
slightly displacing the feed from the reflector focus. Defocusing can cause a degradation 
at one frequency but an improvement at the other Relative to an optimum pattern, 
defocusing will result in a spoiled beam; that is, null filling and change in beam slope. An 
example is shown in Figure 19. Defocusing also affects pattern levels in the rear 
hemisphere 

C. OPTIMUM DESIGNS 

1. Single Frequency (1850 MHz) 

After making all the comparisons mentioned earlier, an optimum configuration 
was identified for the 36 inch parabolic reflector operating at 1850 MHz. The patterns are 
shown in Figure 20, and the configuration parameters shown in Table 5. 
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0 




Figure 19. Radiation Pattern with Reflector Shifted 0.2A. in the Positive Direction To 
Defocus the Feed. 
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9, degrees 



Figure 20. Radiation Pattern of Optimum Single Frequency Design for 36 inch Reflector 
(The Dimensions are given in Table 6). 
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TABLE 5. OPTIMUM GEOMETRY FOR 36 INCH 
REFLECTOR AT LOW FREQUENCY 


DM 


= 


4.1 


DNL 


= 


0.48 


DS 


- 


0.5 


PLN 


= 


0.4 


ZC 


= 


0.3 


DDIP 


= 


-0.25 


FOD 


= 


0.4 


ROT 


= 


0 


RC 


= 


0.5 


AW 


= 


0.001 


ZR 


= 


0 


ZSHIFT 


= 


0.2 


SW 


= 


0.05 









This radiation pattern meets almost all the requirements listed in the Appendix 
A. The calculated total gain is G = 20.7775, which is sufficient, but the far out sidelobes 
are about 2 dB too high. Unfortunately this is a narrow band design and has poor 
performance at the lower frequency (1350 MHz). This design illustrates that achieving 
the specifications at a single frequency is relatively easy. However, broadbanding the 
antenna proves to be much more difficult. 

2. Dual Frequency (1350 and 1850 MHz) 

Starting from the configuration described in Table 4, attempts have been made to 
obtain a broadband design for the frequency range from 1350 to 1850 MHz. Figures 21 
and 22 present the radiation patterns obtained for the 36 inch antenna which achieves 
most of the specifications given in the Appendix A. The gains are marginal at both 
frequencies and the beamwidth at the low frequency in too wide. Table 6 describes the 
configuration of this reflector. 



38 




Figure 21. Radiation Pattern of 36 inch Reflector for Lower Frequency (1350 MHz). 
Dimensions are given in Table 6. 
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0 




Figure 22. Radiation Pattern of 36 inch Reflector for Higher Frequency (1850 MHz). 
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TABLE 6. OPTIMUM GEOMETRY FOR 36 INCH 
REFLECTOR FOR THE LOW FREQUENCY 


DM 


= 


4.1 


DNL 


— 


0.5 


DS 


= 


0.6 


PLN 


= 


0.4 


ZC 


= 


0.3 


DDIP 


- 


-0.25 


FOD 


= 


0.4 


ROT 


= 


0 


RC 


= 


0.6 


AW 


= 


0 


ZR 


= 


0.1 


ZSHIFT 


= 


0.25 


SW 


= 


0.05 









For the 48 inch antenna. Figures 23 and 24 show the radiation patterns for the 
optimum configuration whose dimensions are listed in Table 7. 
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Figure 23. Radiation Pattern of 48 inch Reflector for Lower Frequency (1350 MHz). 
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0 




Figure 24. Radiation Pattern of 48 inch Reflector for Higher Frequency (1850 MHz). 
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TABLE 7. OPTIMUM GEOMETRY FOR 48 INCH 
REFLECTOR FOR HIGH FREQUENCY 


DM 


= 


7.5 


DNL 


= 


0.48 


DS 


= 


0.8 


PLN 


= 


0.4 


ZC 


= 


0.3 


DD1P 


= 


-0.25 


FOD 


= 


0.4 


ROT 


= 


0 


RC 


= 


0.5 


AW 


= 


0.001 


ZR 


= 


0 


ZSHIFT 


= 


0.2 


SW 


= 


0.05 









It is worth noting that there is no feed configuration that can be used for both frequencies 
and simultaneously meet the specifications. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



The objective of this thesis was to find the optimum combination of the feed and 
reflector geometry over the frequency band 1350 to 1850 MHz. The current antenna size 
is 36 x 48 inches. Since the computer code only accepts circular reflector apertures, two 
separate circular antennas were analyzed One was 36 inches and the second was 48 
inches. The proposed feed design, which consists of a fed dipole in the vicinity of a 
parasitic ring and disk, was a good approach to achieve a significant improvement 
compared with the original geometry of the Marine radio AN/MRC-142. Because of its 
simplicity, the cost of the new feed represents a minimum impact on the system upgrade. 

Calculations for both sizes were made separately and, as it is shown in Chapter III, 
the 36 inch reflector could not satisfy the requirements at both frequencies 
simultaneously, but it gave the best result for the high frequency ( 1 850 MHz). The results 
obtained with the 48 inch reflector come very close to satisfying all the requirements, 
although the sidelobe level (SLL) is off by approximately 2 dB at wide angles in the 
H-plane. 

An assumption with regard to the transmission line losses was made, which affect 
the total gain of the antenna In the field calculations, solid surface reflection 
characteristics were assumed. In fact, the surface is a grid, and consequently the 
reflection coefficients for the 0 and <j) components of the electric field will vary with the 
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angle of incidence. This effect will be small if the wire spacing is small compared to X. 
Achieving the required gain was not a problem for the 48 inch design, and small 
adjustments in the antenna feed geometry did not have a significant negative impact on 
the antenna performance. It was also determined that if there was an improvement due to 
an additional parasitic dipole, the increase in complexity did not merit its use. 

To summarize, a low sidelobe reflector antenna has been designed which fulfilled 
most of the requirements for the Marine radio. The dipole disk feed with beam-forming 
ring is electrically and mechanically simple It can be supported by means of a stiff 
coaxial cable aligned along the axis of the paraboloid, and held by a center support which 
offers minimum blockage when compared with the usual strut supports, as was stated by 
Kildal and Skyttemyer [Ref. 3], 
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APPENDIX A. ANTENNA SPECIFICATIONS 



Over the frequency range of 1350-1850 MHz, the antenna must have gains as 
described in Figure A. 




Azimuth Relative to Boresight (degrees) 


Gain (dBi) 


AtO 


A - see note A 


At± 10 


At least 5 dB lower than A 


At ±20 


At least 19.5 dB lower than A 


At ±30 to 180 


At least 27 dB lower than A 



Shall satisfy the relationship 17.5 < A < 21 dB 



Figure A. Gain Requirements. 
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TABLE A. ANTENNA/MAST REFERENCE DATA 


Antennna 


Lightweight, open grid, reflector 


Frequency Range 


1350 to 1850 MHz 


Gain 




Frequency 




1350 MHz 


18 dBi minimum 


1850 MHz 


20 dBi minimum 


Breamwidth 


+7.5° 


Horizontal 




7Vz ±1° from boresight 


3 dB half power 


15 ± 1° from boresight 


15 dB down 


Sidelobes 


10 dB min 


Antennna 




Vertical 




V/2 ± 1°, -2° from boresight 


3 dB half power 


15 ± 1°, -2° from boresight 


15 dB down 


Sidelobes 


10 dB min 


Polarization 


Horizontal or vertical polarization 


Weigth 


17.5 lbs 


Maximum Input Power 


15 W 


Maximum Output Power 


2320 W Effective Radiated Power (ERP) 


Reflector Dimensions 


45 in X 35 in. Depth: 16 in 


Mast 


Telescoping mast, consists of eleven 
carbon-fiber tubular sections 


Height 


Continously adjustable from 5 ft 8 in to 50 ft 


Length (for storage and transit) 


5 ft 8 in 


Maximum Height (deployed) 


50 ft 


Headload 


37 lbs maximum 


Maximum Wind (50 ft. erection) 


60 mph (without ice) 

40 mph (with Vi in ice loading) 


Weight 


35 lbs 


2-Man Deployment Time 


15 minutes 


Mast Accessories Time 


Stowage of mast accessories 




Total accessory bag weight: 44 lbs 
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APPENDIX B. PROGRAM CODE* 



Gain2f 

C PROGRAM mrcint.f 
C 

C FAR FIELD PATTERN INTEGRATION FOR GAIN 
CALCULATIONS. 

C *** general wireybor geometries can be handled *** 

C 

C »» INPUT FILE IS pcurrent - SPECIALIZED FOR 
PROGRAM mrc.f 

C »» HAS A GAIN NORMALIZATION LOOP 
INCORPORATED 

C »» WIRE DOES NOT HAVE TO LIE IN THE Z=0 

PLANE 

C 

DIMENSION 

RH(1000),ZH(1000).AG(100),XG(100),XT(10),AT(10) 

DIMENSION 

Q(30),S(30),nw1(2),nw2(2),pha(l000),xh(l000),yh(1000) 

COMPLEX 

RS(1800),ETW,EPW,U,CNT(1500),ET1 ,EP1,ET2,EP2 
COMPLEX 

CONJG,CMPLX,EXP1 ,EXP2,RW(2000).ET,EP,UC 
DATA pi/3.14159/ 
data nt/2/,xt( 1 ) ,at( 1 )/.5773 503 , 1 ./ 
c data 

nt/4/,xt( 1 ) ,xt( 1 0) ,at( 1 ) ,at( 1 0)/.33998 1 04 , .861 1 363 1 1 5 , 
c * .6521451548,-3478548451/ 

RAD=PI/180. 

BK=2.*PI 
U=(0.,1 .) 

UC=-U/4./PI 

OPEN(1 ,FILE=’outgaus’) 

READ(1 ,*) NPHI 
NT2=NT/2 
DO 1 K=1 ,NT2 
K1=NT-K+1 
AT(K1)=AT(K) 

XT(K1)=XT(K) 

1 XT(K)=-XT(K) 

DO 3 K=1 ,NPHI 
3 READ(1 ,*) XG(K),AG(K) 
OPEN(3,file=’pcurrent’,STATUS=’old') 

READ(3,*) 

DM,AW,ZC,RC,DS,NP1 ,NP2,NBLOCK,NWIRES,FFAC 
READ(3,*) 

NP.MT1 ,MP1 ,N.MT2,NSURF,NROW,MHl, MODES 

write(6,*) ’NP1,NP2,N, MODES, NROW,NWIRES,NP=’ 

wrrte(6,*) NP1 ,NP2,N, MODES, NROW,NWIRES,NP 

READ(3,*) (xh(i),i=1 ,np) 

READ(3,*) (yh(i),i=1 ,np) 

READ(3,*) (zh(i),i=1 ,np) 

READ(3,*) (rh(i),i=1.np) 

READ(3,’) (pha(i).i=1,np) 

READ(3,*) (nw1(i),i=1,nwires) 

READ(3,*) (nw2(i),i=1 ,nwires) 



READ(3,*) (cnt(i),i=1,nrow) 

CLOSE(3) 

NS=NBLOCK*N 

OPEN(8.FILE=’outint') 

WR1TE(8,8000) NP1 .NP2.NBLOCK 
8000 FORMAT(//,10X,’** WIRE/BOR ANTENNA GAIN 
CALCULATION 

7/.10X. ’NUMBER OF BOR POINTS (NP1)=’,i4,/,10x, 
•’NUMBER OF DIPOLE POINTS (NP2)=’,I4,/,10X, 
•'NUMBER OF MODE BLOCKS (2-M+1)=’,l2) 

NDIVT=2 

NDIVP=2 

WR!TE(6,*) ’ENTER NDIVT AND NDIVP' 

READ(5,*) NDIVT, NDIVP 
S1=0. 

S2=180. 

Q1=0. 

Q2=180. 

WRITE(6,*) THETA START AND STOP ANGLES 
(DEGREES)’ 

READ(5,*) SI ,S2 

WRITE(6,*) THETA START AND STOP ANGLES 
(DEGREES)' 

READ(5,*) Q1,Q2 
C 

C DEFINE STANDARD BOR INDICES 
C 

MTl=NPl-2 

MP1=NP1-1 

MT2=NP2-2 

C NWIRES=1 SINCE DUMMY SEGS WILL HAVE ZERO 
CURRENT 
NWIRES=1 
NW1(1)=NP1+1 
NW2(1)=NP 
C 

C GAIN NORMALIZATION LOOP (USE THETA=PHI=0) 

C 

ET1=(0.,0.) 

ET2=(0.,0.) 

ETWs(0.,0.) 

CALL 

PLANEW(NWIRES,NW1 ,NW2,NP1 ,NP2,xh,yh,zh,0.,0..RW) 
DO 610 L=1,MT2 

610 ETW=ETW+RW(L)*CNT (L+NS) 

EXP1=(1.,0.) 

EXP2=CONJG(EXP 1 ) 

DO 611 M=1 ,MHI 
NM=M-1 

CALL PLANES(NM,NP1,0,NT,RH,ZH,XT,AT,0.,RS) 

NTOP1=MODES-NM 

NTOP2=NBLOCK*(NTOP1+1) 

NS2=NTOP1*N 
NS1=NTOP2*N 
DO 612 L=1 ,MT 1 

ET1=ET1+RS(L)*CNT(L+NS1)*EXP1 
IF(NM.EQ.O) GOTO 612 
ET1=ET1+RS(L)*CNT (L+NS2)*EXP2 



* Although this is FORTRAN code, the column conventions are not observed for 
this listing so that space can be conserved. 
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612 CONTINUE 
DO 613 L=1 ,MP1 

ET2=ET2+RS(L+MT 1 )*CNT(L+NS1 +MT1 )*EXP1 
IF(NM.EQ.O) GO TO 613 
ET2=ET2*RS(L+MT 1 )*CNT(L+NS2+MT1 )*EXP2 

613 CONTINUE 
611 CONTINUE 

ET=UC*(ET1 +ET2+ETW) 

EMAX=CABS(ET) 

WRITE(6,*) ’EMAX=’,EMAX 

WRITE(8,8008) NT,NPHI,NDIVP,NDIVT,S1 ,S2,Q1 .02 
8008 FORMAT(//,10X,’NT=',l3,/,10x,’NPHI=’,l4,/,10X, 

* ’NUMBER OF INTEGRATION INTERVALS IN 
PHI=',i3,/,10x, 

* 'NUMBER OF INTEGRATION INTERVALS IN 
THETA=',i3,/,10x, 

* THETA LIMITS: START, STOP=’,F7.2,2X,F7.2,/,1 OX, 

* ’PHI LIMITS: START, STOP=’,F7.2,2X,F7.2) 
c 

c theta integration intervals 
c 

S1=S1*RAD 

S2=S2*RAD 

Q1=Q1*RAD 

Q2=Q2*RAD 

DS=(S2-S1)/FLOAT(NDIVT) 

DO 502 1=1 .NDIVT+1 
S(I)=(I-1)*DS 
WRITE(6,’) ’l,S(l)=M,S(l) 

502 CONTINUE 
c 

c phi integration intervals 
c 

DQ=(Q2-Q1)/FLOAT(NDIVP) 

DO 503 1=1 .NDIVP+1 
Q(I)=(I-1)*DQ 
WRITE(6/) ■|,Q(I)=M,Q(I) 

503 CONTINUE 
SUM=0. 

do 2000 JJ=1 .NDIVP 
P1=DQ/2. 

P2=(Q(JJ+1 )+Q(JJ))/2. 

DO 2000 J=1 ,NPHI 
PHR=p1*xg(j)+p2 
DO 1000 ll=1,NDIVT 
T1=DS/2. 

T2=(S(ll+1)+S(ll))/2. 

DO 1000 1=1 ,NPHI 
THR=T1 *XG(I)+T2 
IF(THR.LT.PI) GO TO 510 
THR=(BK*THR) 

PHR=PHR+PI 
510 CONTINUE 
STRA=SIN(THR)*AG(l)*ag(j) 

ET1=(0.,0.) 

EP1=(0.,0.) 

ET2=(0.,0.) 

EP2=(0.,0.) 

ETW=(0.,0.) 

EPW=(0.,0.) 

C 

C WIRE CONTRIBUTION TO THE RADIATED FIELD 
C 



CALL 

PLANEW(NWIRES,NW1 ,NW2,NPl ,NP2,XH.YHZH,THR,P 
HR,RW) 

DO 210 L=1 ,MT2 
ETW=ETW+RW(L)*CNT(L+NS) 

210 EPW=EPW+RW(L+MT2)*CNT(L+NS) 
c 

c BOR SURFACE CONTRIBUTION TO THE RADIATED 

FIELD 

c 

DO 300 M=1 ,MHI 
NM=M-1 

EXP1=CEXP(CMPLX(0.,FLOAT(NM)’PHR)) 

EXP2=CONJG(EXP1) 

CALL PLANES(NM,NP1 ,0,NT,RH2H,XT,AT,THR,RS) 

NTOP1=MODES*NM 

NTOP2=NBLOCK-(NTOP1 + 1) 

NS2=NTOP1*N 
NS1=NTOP2’N 
DO 250 L=1,MT1 

ET1=ET1+RS(L)*CNT(L+NS1)*EXP1 
EP1=EP1+RS(L+N)*CNT(L+NS1)*EXP1 
IF(NM.EQ.O) GOTO 250 
ET1=ET1+RS(L)*CNT(L+NS2)*EXP2 
EP1=EPTRS(L+N)*CNT(L+NS2)*EXP2 
250 CONTINUE 
DO 260 L=1 ,MP1 

ET2=ET2+RS(L+MT1 )*CNT(L+NS1+MT1)*EXP1 
EP2=EP2+RS(L+MT1+N)*CNT(L+NS1+MT1)*EXP1 
IF(NM.EQ.O) GOTO 260 
ET2=ET2-RS(L+MT 1 )*CNT(L+NS2+MT1 )*EXP2 
EP2=EP2+RS(L+MT1+N)*CNT(L+NS2+MT1)*EXP2 
260 CONTINUE 
300 CONTINUE 
ET=UC*(ET1+ET2+ETW) 

EP=UC*(EP1+EP2+EPW) 
SUM=SUM+(CABS(ET)**2+CABS(EP)**2)’STRA 
1000 CONTINUE 
2000 CONTINUE 
PRAD=T1’P1*SUM 
GAIN=4.*Pl*EMAX**2/PRAD 
GDB=10.*alog1 0(GAIN) 

WRITE(8,8010) PRAD,EMAX,GAIN,GDB 
8010 FORMAT(//,1 OX, TOTAL RADIATED 
POWER=’,E15.8,/,10X, 

* ’MAXIMUM FIELD VALUE (V/M)=’,E15.8,/,10X, 

* ’GAIN=’,E1 5.8,/, 1 0x,’(dB)=’,f8.2) 

WRITE(6,*) TRAD, GAIN, dB=\PRAD, GAIN, GDB 

9999 STOP 
END 

SUBROUTINE 

PLANES(M1,NP1,NP2,NT,RH,ZH,XT,AT,THR,R) 

C 

C PLANE WAVE EXCITATION VECTOR FOR BOR 
ELEMENTS. 

C NO CHANGE FROM HARRINGTON’S ORIGINAL 
PROGRAM OTHER THAN 
C *• ONLY ONE MODE PER CALL ” 

C 

COMPLEX 

R(1800),U t U1,UA,UB,FA(50) f FB(50), F2A.F2B.F1 A, FIB 
COMPLEX U2,U3,U4,U5,CMPLX 
DIMENSION RH( 1 000) ,ZH( 1 000) ,XT( 1 0) ,AT( 1 0) 
DIMENSION R2( 1 0) ,Z2( 1 0) ,BJ(5000) 

U=(0.,1.) 
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U1=3.141593*U**M1 

NF=1 

M2=M1 

NP=NP1+NP2 

MP=NP-1 

MT=MP-1 

N=MT+MP 

N2=2*N 

CC=COS(THR) 

SS=SIN(THR) 

M3=M1 + 1 

M4=M2+3 

IF(MI.EQ.O) M3=2 

M5=M1+2 

M6=M2+2 

DO 12 IP=1 ,MP 

K2=IP 

l=IP+1 

DR=.5*(RH(I)-RH(IP)) 

DZ=.5*(ZH(I)-ZH(IP)) 

D1=SQRT(DR*DR+DZ*DZ) 

R1=.25*(RH(I)+RH(!P)) 

IF(ABS(R1).LT.1 .E-5) R1=1. 

Z1=.5*(ZH(I)+ZH(IP)) 

DR=.5*DR 
D2=DR/R1 
DO 13 L=1 ,NT 
R2(L)=R1+DR*XT(L) 
Z2(L)=Z1+DZ*XT(L) 

13 CONTINUE 
D3=DR*CC 
D4=*DZ*SS 
D5=DTCC 
DO 23 M=M3,M4 
FA(M)=0. 

FB(M)=0. 

23 CONTINUE 
DO 15 L=1 ,NT 
X=SS*R2(L) 

IF(X.GT..5E-7) GOTO 19 
DO 20 M=M3,M4 
BJ(M)=0. 

20 CONTINUE 
BJ(2)=1 . 

S=1. 

GOTO 18 
19 M=2.8*X+14.*2./X 
IF(X.LT..5) M=11.8+ALOG10(X) 
IF(M.LT.M4) M=M4 
BJ(M)=0. 

JM=M*1 
BJ(JM)=1 . 

DO 16 J=4,M 
J2=JM 
JM=JM*1 
J1=JM-1 

BJ(JM)=J1/X*BJ(J2)-BJ(JM+2) 

16 CONTINUE 
S=0. 

IF(M.LE.4) GOTO 24 
DO 17 J=4,M,2 
S=S+BJ(J) 

17 CONTINUE 

24 S=BJ(2)+2.*S 

18 ARG=Z2(L)*CC 



UA=AT(L)/S’CMPLX(COS(ARG),SIN(ARG)) 
UB=XT{L)*UA 
DO 25 M=M3,M4 
FA(M)=BJ(M)*UA+FA(M) 

FB(M)=BJ(M)*UB+FB(M) 

25 CONTINUE 
15 CONTINUE 

IF(MI.NE.O) GOTO 26 
FA(1)=* FA(3) 

FB(1)=- FB(3) 

26 UA=U1 

DO 27 M=M5,M6 

M7=M*1 

M8=M+1 

F2A=UA*(FA(M8)+FA(M7)) 

F2B=UA*(FB(M8)+FB(M7)) 

UB=U*UA 

FI A=UB*(FA(M8)-FA(M7)) 

F1B=UB*(FB(M8)-FB(M7)) 

U4=D4*UA 

U2=D3*F1 A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT)=*D5*(F2A+D2*F2B) 

R(K5+MT)=D1*(F1 A+D2*F1B) 

IF(IP.EQ.I) GO TO 21 
R(K1)=R(K1)+U2-U3 
R(K4)=R(K4)+U4-U5 
IF(IP.EQ.MP) GO TO 22 

21 R(K2)=U2+U3 
R(K5)=U4+U5 

22 K2=K2+N2 
UA=UB 

27 CONTINUE 
12 CONTINUE 

RETURN 

END 

SUBROUTINE 

PLANEW(NWIRES,NW1 ,NW2,NP1 ,NP2,XH,YH,ZH, 

* THR,PHR,R) 

C 

C PLANE WAVE EXCITATION VECTOR ELEMENTS FOR 
WIRE AND 

C ATTACHMENT SEGMENT. INCIDENCE DIRECTION IS 
(THR.PHR). 

C * ONLY ONE INCIDENCE ANGLE PER CALL 

C * WIRE DOES NOT HAVE TO LIE IN THE Z=0 PLANE 

C 

COMPLEX 

U0 1 AA.BB,C,R(400),CEXP,EXP,FI1 1 FI2,SI,DI,CMPLX 

DIMENSION 

ZH(400),NS(2) 1 NW1(2),NW2(2),XH(400),YH(400) 
MP2=NP2-1 
MT2=NP2*2 
MT22=2*MT2 
DO 5 L=1,NWIRES 
5 NS(L)=NW2(L)-NW1 (L)+1 
U0=(0.,0.) 

CC=COS(THR) 

SS=SIN(THR) 

CP=COS(PHR) 
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SP=S1N(PHR) 

UP=SS*CP 

VP=SS*SP 

WP=CC 

DO 12 IP=1,MP2 

II=1P+NP1 

1 = 11+1 

ZS=.5*(ZH(I)+ZH(II)) 

XS=.5*(XH(l)+XH(ll)) 

YS=.5*(YH(I)+YH(II)) 

DX=XH(l)-XH(ll) 

DY=YH(l)-YH(lt) 

D1=SQRT(DX*'2+DY**2) 

SU=DY/D1 

CU=DX/D1 

C FOR WIRES IN PERPENDICULAR TO Z SIN(V)=1 AND 
COS(V)=0 
SV=1.0 
CV=0.0 
C 

C WIRE SEGMENT CALCULATIONS 
C 

A=UP*CU+VP*SU 

B=UP*XS+VP*YS+WP*ZS 

C=CMPLX(0.,A) 

EXP=CEXP(CMPLX(0.,B)) 

AA=CC*(CU*CP+SU*SP) 

BB=SU*CP-SP*CU 
PSI=D1 *A'2. 

IF(PSI.NE.O.) GOTO 60 
SINC=1 . 

GO TO 61 

60 SINC=SIN(PSI)/PSl 

61 COSP=COS(PSI) 

Fl1=SlNC*D1*EXP/2. 

FI2=(0.,0.) 

IF(ABS(A).LT. 1 .E-4) GO TO 62 
CSP=COSP-SlNC 
IF(ABS(CSP).LT.1 .E-4) GO TO 62 
FI2=EXP/C’CSP 

62 CONTINUE 
SI=F11+FI2 
DI=FI1-FI2 

C 

C R-WIRE-THETA 
C 

IF(IP EQ.MP2) GOTO 10 
R(lP)=AA*Sl 
R(1P+MT2)=BB*SI 
10 CONTINUE 
C 

C R-WIRE-PHl 
C 

14 IF(IP.EQ.I) GOTO 12 
R(IP-1)=R(IP-1 )+AA*Dl 
R(lP-1+MT2)=R(lP-1+MT2)+BB*Dl 
12 CONTINUE 
RETURN 
END 

Reflector 

C PROGRAM mrc.f (Marine Radio Antenna Program) 

C 

c Modified from ORDONEZ.f -- hyperboloid becomes a disc 



C Paraboloie WITH A PARASITIC RING DIPOLE FEED. 
Fed dipole is at z=0. 

C This program uses a parasitic dipole and two h-plane 
C parasitic dipoles and the subreflector. This system also 
contains 
C a cavity ring, 

C 

C ICALC=0 CALC CURRENTS AND FIELD 
C =1 CALC FIELDS ONLY (READ CURRENTS) 

C IMP=0 PEREECTLY CONDUCTING SURFACES 
C =1 SOME NONZERO SURFACE RESISTANCE 
C ICWRT=0 WRITE CURRENT COEFFICIENCTS TO 
FILE pcurrent 

C IRES=0 READ RESISTIVE CORRECTIONS 
C 

COMPLEX 

EP,ET,Z(500000),RS(1000),RW(400),B(700),C(700),U,UC 
COMPLEX ZL0(400),ZL(1500),EC,EX,ZSS(50000) 
COMPLEX EXP1 ,EXP2,ZSW(50000),ZWW( 10000) 
COMPLEX 

CEXP,CONJG,CMPLX,ET1 ,ET2,EP1 ,EP2,ETW,EPW 
DIMENSION 

RH(400),ZH(400),XT(4),AT(4),IPS(700),NW1 (4),NW2(4) 
DIMENSION 

A(400),X(400).EXP(800) t ANG(800),ECP(800) 

DIMENSION ECV(800),EXV(800),PHC(800),PHX(800) 
DIMENSION 

XX(41),YY(41),XH(400).YH(400),PHA(400) 

DATA 

PI, START, STOP, DT/3.14159,0., 180., 1./.ICWRT/0/, IRES/1/ 
DATA lCALC/0/,IPRINT/0/,IMP/0/,lTEST/1/ 

DATA nt/2/,xt( 1 ) ,at( 1 )/.5773503,1 ./ 
c data nt/4/,xt(1),xt(2),at(1),at(2)/. 339981 04, .861 1363115, 
c * .6521451548, .3478548451/ 

OPEN(1,FlLE='outgaus’,status= , old') 

OPEN(30,FlLE=’rh.m’) 

OPEN(31 ,FILE=’zh.m’) 

READ(1 ,*) NPHI 
DO 3 K=1 ,NPHl 
READ(1,*) X(K),A(K) 

3 CONTINUE 
RAD=PI/180. 

BK=2.*PI 

U=(0.,1.) 

U0=(0.,0.) 

UC=-U/4./PI 
NT2=NT/2 
DO 1 K=1 ,NT2 
K1=NT-K+1 
AT(K1)=AT(K) 

XT(K1)=XT(K) 

1 XT(K)=-XT(K) 

C 

C NP10 = NUMBER OF PARABOLOID POINTS 
C NP1 1 = NUMBER OF DISC POINTS 
C NP12 = NUMBER OF RING POINTS 
C 

MODES=1 

MHI=MODES+1 

NBLOCK=2*MODES+1 

NWIRES=1 

WRITE(6,*) ’ENTER START. STOP, DT’ 

READ(5,*) START, STOP, DT 

NP10=61 

NP1 1=31 
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NP12=3 

WRITE(6,*) ’ENTER PARA (NP10), DISC (NP1 1) AND 
RING (NP12) PTS' 

READ(5,*) NP10.NP1 1,NP12 
NP1=NP10+NP1 1+NP12 
C 

C DM=MAIN REFLECTOR DIAMETER 
C DS=DISC (SUBREFLECTOR) DIAMETER 
C ZC=DISC LOCATION ALONG Z AXIS 
C FOD=f/DM RATIO 
C RC=RADIUS OF THE RING 
C SW=RING (STRIP) WIDTH 
C 

DM=5.0 

R0=DM/2. 

DS=1 .1 
ZC=.25 
FOD=.3 
RC=0.5 
SW=.05 

WRITE(6,*) ’ENTER DM, DS, ZC (>0), FOD, RC, ZR 
AND SW’ 

READ(5,*) DM,DSZC,FOD,RC,ZR,SW 
FM=FOD*DM 

WRITE(6,*) ’ENTER ZSHIFT FOR PARABOLOID’ 
READ(5,*) ZSHIFT 

GENERATE PARABOLOID CONTOUR 

PHIV=2.’ATAN(1 ./4./FOD) 

FNM=FLOAT(NP10-1) 

DO 50 l=1,NP10 

TH=FLOAT(l-1)*PHIV/FNM 

RM=2.*FM/(1.+COS(TH)) 

XH(l)=0. 

YH(l)=0. 

ZH(l)=-RM*COS(TH)’BK+ZSHIFT*BK 

RH(I)=RM*SIN(TH) , BK 

PHA(ll)=0. 

50 CONTINUE 

DISC LOCATED A Z=ZC 

DSK=DS/2./FLOAT(NP1 1-1 ) 

DO 5 1=1 ,NP1 1 

ll=l+NP10 

XH(ll)=0. 

YH(ll)=0. 

ZH(II)=ZC*BK 

RH(ll)=FLOAT(l-1)*DSK*BK 

PHA(ll)=0. 

CONTINUE 

NEXT THE SIDEWALLS 
C 

DSW=SW/FLOAT(NP12-1) 

ZRBK=ZR*BK 
IF(NP12.NE.O) THEN 
DO 7 1=1 ,NP12 
ll=NP1 1+NP10+I 
XH(ll)=0. 

YH(ll)=0. 

RH(II)=RC*BK 

ZH(ll)=FLOAT(l-1 )*DSW*BK+ZRBK 
7 PHA(ll)=0. 



ENDIF 

NP1=NP10+NP1 1+NP12 
C 

C DIPOLE ANTENNA POINTS IN WAVELENGTHS: AW= 
RADIUS; NP2= NUMBER 

C OF POINTS (=NP21+NP22; NP21 MUST BE ODD); 
DLGTH= LENGTH; 

C PLGTH IS PARASITIC ELEMENT LENGTH. 

C 

AW=.001 

DLGTH=.45 

PLGTH=.44 

WRITE(6 t *) 'ENTER DLN, PLN, DDIP ROT AND AW’ 
WRITE(6,’) ’(DDIP<0 IS TOWARD PARABOLOID)’ 
READ(5,*) DLGTH, PLGTH, DDIP, ROT,AW 
CENT=DLGTH/2. 

AK=AW*BK 
NP21 = 1 1 
NP22=0 

WRrTE(6.*) ’ENTER NP21 AND NP22 (ONLY 1 DIP IF 
NP22=0)’ 

READ(5,*) NP21,NP22 
NW1(1)=NP1 + 1 
NW2(1)=NP1+NP21 
NP2=NP21 +NP22 
DEL=DLGTH/FLOAT(NP21 - 1 ) 

DO 10 1=1 ,NP21 

ll=NP1+l 

YH(ll)=0. 

ZH(ll)=0. 

XH(ll)=BK*(FLOAT(l-1)*DEL-CENT) 

RH(II)=SQRT(XH(II)*’2+YH(II) ,, 2) 

PHA(ll)=ATAN2(YH(ll),XH(ll)+1.E-6) 

10 CONTINUE 
IF(NP22.NE.O) THEN 
NWIRES=2 

NW1 (2)=NP 1 +NP21 + 1 

NW2(2)=NP1+NP21+NP22 

DEL=PLGTH/FLOAT(NP22-1) 

CENT=PLGTH/2. 

DO 11 1=1 ,NP22 
ll=NP1+NP21+l 
ZH(II)=DDIP*BK 

XH(ll)=BK*(FLOAT(l-1) , DEL-CENT)*COS(ROT*RAD) 
YH(ll)=BK*(FLOAT(l-1 )*DEL-CENT)*SIN(ROT*RAD) 
RH(II)=SQRT(XH(II)**2+YH(II)**2) 
PHA(ll)=ATAN2(YH(ll),XH(ll)+1.E-6) 

11 CONTINUE 
ENDIF 

C NP23=NUMBER OF H-PLANE DIPOLE POINTS 
C HLGTH=H-PLANE DIPOLE LENGTH 
C DHP=H-PLANE DIPOLE DISTANCE FROM FEED 
C ZHP=H-PLANE DIPOLE Z COORDINATE 
C NP24=NUMBER OF POINTS FOR THE OTHER 
H-PLANE DIPOLE 
NP23=7 
NP24=7 
HLGTH=.2 
DHP=.2 
ZHP=.001 

WRrTE(6,*) 'ENTER NUMBER OF POINTS PER 
H-PLANE DIPOLE' 

WRfTE(6,*) ’(0 FOR NO H-PLANE PARASITIC 
DIPOLES)’ 

READ(5,*) NP23 
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IF(NP23.NE.O) THEN 

WRITE(6,*) 'H-PLANE DIPOLE LENGTH, DISTANCE 
AND LOCATION’ 

READ(5,‘) HLGTH.DHP^HP 

NWIRES=NWIRES+2 

NP24=NP23 

NP2=NP21 +NP22+NP23+NP24 
NW1 (3)=NP1+NP21+NP22+1 
NW2(3)=NP1+NP21+NP22+NP23 
NW1 (4)=NP1 +NP21 +NP22+NP23+1 
NW2(4)=NP1+NP21+NP22+NP23+NP24 
DEL=HLGTH/FLOAT(NP23-1) 

CENT=HLGTH/2. 

DO 12 1=1 ,NP23 
I1=NP1+NP21+NP22+1 
2H(I1)=2HP*BK 
YH(I1)=DHP*BK 

XH(l1)=BK*(FLOAT(l-1)*DEL-CENT) 

RH(I1)=SQRT(XH(I1)**2+YH(I1)**2) 

PHA(l1)=ATAN2(YH(l1),XH(l1)+1.E-6) 

12=11 +NP23 

2H(I2)=2HP*BK 

YH(I2)=-DHP*BK 

XH(l2)=BK*(FLOAT(l-1)*DEL-CENT) 
RH(I2)=SQRT(XH(I2)**2+YH(I2)**2) 
PHA(I2)=ATAN2(YH(I2),XH(I2)+1 E-6) 

12 CONTINUE 
ENDIF 
CLOSE(I) 

NP=NP1+NP2 

WRITE(6,*) 'ENTER FREQUENCY FACTOR’ 

READ(5,*) FFAC 
OPEN(8,file='outmrc’) 

WR ITE(8,8000) 

DM,D$,FOD,PHIV/RAD,NP10,NP1 1 ,RC,SW, 

NP12,2SHIFT,DLGTH,PLGTH,DDIP,NP2,AW,MODES,NT, 

NPHI.FFAC 

IF(NP23.NE.O) WRITE(8,8001) HLGTH,DHP,2HP,NP23 

8000 FORMAT(//,5X,’*** MARINE RADIO REFLECTOR 
ANTENNA 5X, 

‘’REFLECTOR PARAMETERS (LENGTHS IN 
WAVELENGTHS):’,/, 5X, 

•’MAIN REFL DIA=’,F8.3,/,5X,'DISC 
DIA=’,F8.3,/,5X,'F/D=',F8.4, 

• /,5X,’PHIV (DEG)=’,F8.2,/,5X,’MAIN REFL POINTS 
(NP1 1)=’,I4,/, 

* 5X/SUB REFL POINTS (NP12)=M4,/,5X,/,5X, 

• ’RING RADIUS=’.F8.3,/,5X, ’SIDEWALL 
LENGTHS, F8.3./.5X, 

• ’NUMBER OF RING POINTS (NP12)=’,I4,//,5X, 

* 'PARABOLOID SHIFT DISTANCE, 

2SHIFT=’,F8.3,/,5X, 

• 'DIPOLE PARAMETERS:’, /,5X, ’FED DIPOLE 
LENGTH=\F8.3,/,5X, 

* ’PARASITIC DIPOLE 
LENGTH=’.F8.3,/,5X,’SPACING=’,F8.3,/,5X, 

• ’NUMBER OF DIPOLE POINTS =’,14, /,5X, ’RADIUS OF 
DIPOLE=’. 

* F8. 4, //,5X, ’NUMBER OF A2IMUTHAL 
M0DES=’,I3,/,5X,’NT=',I3, 

* /,5X,’NPHI=’,I3,/,5X,’FREQ. SCALE FACTOR=’.F8.4) 

8001 FORMAT(/,5X, 'H-PLANE DIPOLE 
LENGTH=’,F8.3,/,5X, 



* ’H-PLANE DIPOLE DISTANCE FROM 
AXIS=’,F8.3,/,5X, 

* 2 POSITION OF H-FIELD 
DIPOLE=’,F8.3,/,5X, ’NUMBER OF H-FIELD 

• DIPOLE POINTS=’,l3) 

IF(ISEG.EQ.O) WRITE(8,1300) 

1300 

FORMAT(/, 5X, ’INDEX’, 8X, 2(1)’, 10X,’RHO(l)',8x,'2SURF’) 
C 

C FILE difdat CONTAINS RESISTIVE CORRECTIONS IF 

IRES.NE.O 

C 

IF((IMP.EQ.0). OR. (IRES.NE.O)) GO TO 70 
OPEN(9,FILE='difdat’) 

READ(9,*) ID 
WRITE(6,*) ’ID=’,ID 
DO 75 1=1, ID 

READ(9,*) ll,RI.YY(l),DUM 
IF(YY(l).LT.0.) YY(l)=0. 

75 CONTINUE 

XXI=1 ./FLOAT(ID-I) 

DO 15 11=1, ID 
XX(ll)=XXI*FLOAT(ll-1) 

15 CONTINUE 

WRITE(6,*) ’SURFACE IMPEDANCE VALUES READ' 
70 CONTINUE 
DO 52 1=1, NP1 

IF(ABS(2H(I)).LT..001) 2H(I)=0. 
IF(ABS(RH(l)).LT..001) RH(l)=0. 

2L0(I)=(0.,0.) 

XH(I)=XH(I)‘FFAC 
YH(I)=YH(I)*FFAC 
2H(I)=2H(I)*FFAC 
RH(I)=RH(I)*FFAC 
2HB=2H(I)/BK 
RHB=RH(I)/BK 
WRITE(30,5019) RHB 
WRITE(31 ,5019) 2HB 
c if(rhb.ge.10.) zl0(i)=(.2,0.) 

C 

C NON2ERO SURFACE IMPEDANCES IF IMP > 0 
C 

IF((IMP.EQ.0). OR. (IRES.NE.O)) GO TO 51 
IF(I.GT.NPIO) GO TO 51 
C 

C USE THE EXACT RESISTIVE CORRECTIONS - 
INTERPOLATE BETWEEN 
C DATA POINTS. 

C 

UU=RHB/R0 

CALL INTERP(UU,VV,XX.YY,ID) 

2L0(I)=VV 

51 CONTINUE 

WRITE(8,8004) l,2HB,RHB,2L0(l) 

52 CONTINUE 

8004 FORMAT(6X,l4,4X,F8.3,8X,F8.3,6X .2F8.4) 
WRITE(8,1310) 

1310 FORMAT(/,5X, 'DIPOLE COORDINATES:’,//, 

• 5X, ’INDEX’, 8X.2(l)’,10X,’RHO(l)M0X, ’PHI PLN’) 

DO 53 l=NP1+1 ,NP 
IF(ABS(2H(I)).LT..001) 2H(I)=0. 
IF(ABS(RH(I)).LT..001) RH(l)=0. 

XH(I)=XH(I)*FFAC 

YH(I)=YH(I)*FFAC 

2H(I)=2H(I)*FFAC 
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RH(I)=RH(I)*FFAC 
PHB=PHA(I)/RAD 
ZHB=ZH(I)/BK 
RHB=RH(I)/BK 
WRlTE(30 f 5019) RHB 
WRITE(31 ,5019) ZHB 
53 WRITE(8,8005) l,ZHB,RHB,PHB 
8005 FORMAT(6X,l4,4X,F8.3,8X,F8.3,9X,F8.2) 

799 CONTINUE 

IF(ITEST.EQ.O) GO TO 9998 
wrlte(6,*) 'geometry defined’ 

C 

C DEFINE DIMENSIONS OF THE IMPEDANCE MATRIX 

BLOCKS 

C 

MT 1=NP1-2 

MP1=NP1-1 

N=MP1+MT 1 

MT2=NP2-2 

NSURF=NBLOCK*N 

NSW=NSURF+MT2 

NROW=NSW 

WRITE(6,*) 

’NP1 ,NP2,NSURF,NROW=’,NP1 ,NP2,NSURF,NROW 
WRITE(6,*) ’MT2 SHOULD BE ODD - MT2=’,MT2 
IF(ITEST.EQ.O) GO TO 9998 
IF(ICALC.EQ.O) THEN 
C 

C COMPUTE IMPEDANCE MATRIX ELEMENTS 
C 

CALL 

ZMATWW(NWIRES,NW1 ,NW2,XH,YH,RH,ZH,NT,XT,AT,A 
K,ZWW) 

CALL ZASMB0(NP1 ,NP2.MODES,Z,ZWW) 
write(6,*) ’wire impedance computed' 

IF(IMP.NE.O) THEN 

CALL ZLOAD(NP1,RH,ZH,ZLO,ZL) 

WRITE(6,*) ’RESISTIVE SURFACE IMPEDANCE 
MATRIX COMPUTED’ 

ENDIF 

DO 400 M=1 ,MHI 

NM=M-1 

CALL 

ZMATSS(NM,NM,NP10,NP1 1 ,NP12,NPHI,NT,RH,ZH,X,A,X 
T.AT.ZSS) 

IF(IMP.NE.O) CALL ZTOT(MT1 ,MP1 ,ZLZSS) 

CALL 

ZMATSW(NWIRES,NW1 ,NW2,NP10,NP1 1 ,NP12,XH,YH,R 
H,ZH, 

* NT,XT,AT,NPHI,X,A,NM,AK,ZSW) 

CALL ZASMBN(NP1 ,NP2, MODES, NM.Z.ZSS2SW) 

IF(NM.EQ.O) GO TO 400 

NMN=-NM 

CALL 

ZM ATSW(N WIRES, NW1 ,NW2,NP1 0,NP1 1 ,NP12,XH,YH,R 
H,ZH, 

* NT,XT,AT,NPHI.X,A,NMN,AK,ZSW) 

CALL ZASMBN(NP1 ,NP2, MODES, NMN2.ZSS2SW) 
400 CONTINUE 
write(6,*) ’z filled’ 

CALL DECOMP(NROW,IPS,Z) 
wrrte(6,*) ’z decomposed’ 

C 

C EXCITATION ELEMENTS (NONZERO FOR THE 
DIPOLE ONLY) 



C 

DO 40 1=1 ,NROW 
40 B(l)=(0.,0.) 

I0=NSURF+(NP21 - 1 )/2 
B(I0)=(1 .,0.) 

CALL SOLVE(NROW,IPS,Z,B,C) 
write(6,*) ’current coefficients computed’ 
c do 665 it=1 ,nrow 
c665 write(6,*) ii,c(ii) 

IF((ICWRT.EQ.0).AND.(ICALC.EQ.0)) THEN 
C 

C WRITE CURRENTS ON DISC FOR PATTERN 
INTEGRATION 
c (program ’cavarchint.f’) 

C 

OPEN(3,file=’pcurrent’) 

WRITE(3,*) 

DM,AWZC,RC,DS,NP1,NP2,NBLOCK,NWIRES,FFAC 

WRITE(3,*) 

NP,MT 1 ,MP1 ,N,MT2,NSURF,NROW,MHI, MODES 
WRITE(3,*) (XH(I),I=1,NP) 

WRITE(3,*) (YH(I),I=1 ,NP) 

WRITE(3,*) (ZH(I),I=1,NP) 

WRITE(3,*) (RH(I),I=1 ,NP) 

WRITE(3,*) (PHA(I),I=1 ,NP) 

WRITE(3,*) (NW1(I),I=1 .NWIRES) 

WRITE(3,*) (NW2(I),I=1 ,NWIRES) 

WRITE(3,*) (C(l),l=1 ,NROW) 

ENDIF 

ENDIF 

C 

C IF ICALC.NE.O THEN READ CURRENT COEFFICIENTS 
FROM DISK FILE 
C 

IF(ICALC.NE.O) THEN 
OPEN(3,file=’pcurrent’,STATUS=’OLD’) 

READ(3,*) 

DMX,AWX,ZCX,RCX.DSX,NP1X,NP2X,NBLX,NWX 

READ(3,*) 

NPX.MT 1X,MP1X,NX,MT2X,NSX,NRX,MHIX,MODESX 
READ(3,*) (xh(i),i=1,np) 

READ(3,*) (yh(i),i=1,np) 

READ(3,*) (zh(i),i=1,np) 

READ(3,*) (rh(i),i=1 ,np) 

READ(3,*) (pha(i),i=1,np) 

READ(3,*) (nw1(i),i=1,nwires) 

READ(3,*) (nw2(i),i=1 ,nwires) 

READ(3,*) (c(i),i=1,nrow) 

CLOSE(3) 

wrrte(6,*) ’data read from pcurrent’ 
do 909 »=1 ,np 

909 wrrte(6,*) ’i,rh,zh=’,i,rh(i)/bk,zh(i)/bk 
ENDIF 

IT=INT((STOP-START)/DT)+1 

C 

C RECEIVER PHI CUTS: DO 0 AND 90 (E- AND H- 
PLANES) 

C 

DO 501 IP=1,2 

ECX=1.e-10 

PHR0=0. 

IF(IP.EQ.2) PHR0=RAD*90. 

DO 500 1=1, IT 

THETA=FLOAT(M)*DT+START 

THR=THETA*RAD 
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PHR=PHRO 

IF(THETA.LE.180.) GO TO 99 
THR=(360.*THETA)*RAD 
PHR=PHRO+PI 
99 CONTINUE 

c wrrte(6,*) ’theta.phis’, theta, phr/rad 
ET 1 =U0 
ET2=U0 
EP1=U0 
EP2=U0 
ETW=U0 
EPW=U0 

DIPOLE FIELD 

CALL 

PLANEW(NWIRES,NW1 ,NW2,NP1 ,NP2,XH,YHZH,THR,P 
HR,RW) 

DO 210 L=1 ,MT2 
ETW=ETW+RW(L)’C(L+NSURF) 

210 EPW=EPW+RW(L+MT2)*C(L+NSURF) 

C 

C REFLECTOR AND CAVITY FIELD CONTRIBUTIONS 
C 

DO 300 M=1,MHI 
NM=M-1 

EXP1=CEXP(CMPLX(0.,FLOAT(NM)*PHR)) 
EXP2=CONJG(EXP1 ) 

CALL 

PLANES(NM,NP10,NP11,NP12,NT,RH,2H I XT 1 AT,THR,RS) 
NTOPl=MODES-NM 
NTOP2=NBLOCK-(NTOP1 + 1 ) 

NS2=NTOP1 *N 
NS1=NTOP2’N 
DO 250 L=1 ,MT1 
ET1=ET1+RS(L)*C(L+NS1 )’EXP1 
EP1=EP1 +RS(L+N)’C(L+NS1)*EXP1 
IF(NM.EQ.O) GO TO 250 
ET1=ET1+RS(L)*C(L+NS2)*EXP2 
EP1=EP1-RS(L+N)*C(L+NS2)*EXP2 
250 CONTINUE 
DO 260 L=1 ,MP1 

ET2=ET2+RS(L+MT1)*C(L+NS1+MT1)*EXP1 
EP2=EP2+RS(L+MT 1+N)*C(L+NS1+MT 1 )*EXP1 
IF(NM.EQ.O) GOTO 260 
ET2=ET2-RS(L+MT1)*C(L+NS2+MT 1 )*EXP2 
EP2=EP2+RS(L+MT 1 +N)*C(L+NS2+MT1 )’EXP2 
260 CONTINUE 
300 CONTINUE 

ET=UC*(ET1+ET2+ETW) 

EP=UC*(EP1+EP2+EPW) 

EC=ET 

EX=EP 

ECV(I)=CABS(EC) 

EXV(I)=CABS(EX) 

ECR=REAL(EC) 

ECI=AIMAG(EC) 

EXR=REAL(EX) 

EXI=AIMAG(EX) 

PHC(l)=ATAN2(ECI,ECR+1.e-10)/RAD 

PHX(l)=ATAN2(EXI,EXR+1.e-10)/RAD 

ANG(I)=THETA 

ECX=AMAX1(ECX,ECV(I),EXV(I)) 

500 CONTINUE 

WRITE(8,103) PHR'RAD.ECX 



103 FORMAT(/,5X, PHI OF RECEIVER 
(DEG)=’,F8.2,/,5X, 

• ’MAXIMUM FIELD VALUE (V/M)=’,E15.5) 

DO 600 1=1, IT 

ECP(I)=(ECV(I)/ECX)**2 

EXP(I)=(EXV(I)/ECX)”2 

ECP(I)=AMAX1(ECP(I),.00001) 

EXP(I)=AMAX1(EXP(I),. 00001) 

ECP(I)=10.’ALOG10(ECP(I)) 

EXP(I)=10.*ALOG10(EXP(I)) 

600 CONTINUE 

IF(IPRINT.NE.O) GO TO 310 
WRITE(8,5015) 

5015 

FORMAT(///, 7X, ’ANGLE’, 15X,’CO-POL\25X t ’X-POL7,7X, 
*’(DEG) , ,4X f 2(’(VOLTS)\4X,’(DEG)\3X, , (DB-REL)\4X)) 
DO 9000 L=1 ,(T 
WRITE(8,5016) 

ANG(L),ECV(L),PHC(L),ECP(L),EXV(L),PHX(L) 

*,EXP(L) 

5016 FORMAT(5X,F6.1 ,3X,2(F8.4,3X,F7.1 ,3X,F7.2,3X)) 
9000 CONTINUE 

310 CONTINUE 
C 

C WRITE MATLAB FILES ( m FILES) 

C 

IF(IP.EQ.I) THEN 
OPEN(2,file='ang.m’) 

OPEN(3,file=’cpole.m’) 

OPEN(4 1 file=’xpole.m’) 

DO 9097 1=1, IT 
WRITE(2,5019) ANG(I) 

WRITE(3,5019) ECP(I) 

9097 WRITE(4,5019) EXP(I) 

CLOSE(2) 

CLOSE(3) 

CLOSE(4) 

ENDIF 

IF(IP.EQ.2) THEN 
OPEN(3,file=’cpolh.m’) 

OPEN(4,file=’xpolh.m’) 

DO 9098 1=1, IT 
WRITE(3,5019) ECP(I) 

9098 WRITE(4,5019) EXP(I) 

CLOSE(3) 

CLOSE(4) 

ENDIF 

501 CONTINUE 
5019 FORMAT(F8.3) 

9998 STOP 
END 

SUBROUTINE SOLVE(N,IPS,UL,B,X) 

COMPLEX UL(500000),B(700),X(700),SUM 
DIMENSION IPS(700) 

NP1=N+1 

IP=IPS(1) 

X(1)=B(IP) 

DO 2 1=2, N 
IP=IPS(I) 

IPB=IP 
IM 1=1-1 
SUM=(0.,0.) 

DO 1 J=1 ,IM1 
SUM=SUM+UL(IP)*X(J) 

1 IP=IP+N 
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2 X(I)=B(IPB)-SUM 
K2=N*(N-1 ) 

IP=IPS(N)+K2 

X(N)=X(N)/UL(IP) 

DO 4 IBACK=2,N 

l=NP1-IBACK 

K2=K2-N 

IPI=IPS(I)+K2 

IP1=I+1 

SUM=(0.,0.) 

IP=IPI 

DO 3 J=IP1,N 
IP=IP+N 

3 SUM=SUM+UL(IP)*X(J) 

4X(I)=(X(I)-SUM)/UL(IPI) 

RETURN 

END 

SUBROUTINE DECOMP(N.IPS.UL) 

COMPLEX UL(500000), PIVOT, EM 
DIMENSION SCL(700),IPS(700) 

DO 5 1=1, N 

IPS(I)=1 

RN=0. 

J1=l 

DO 2 J=1,N 

ULM=ABS(REAL(UL(J1)))+ABS(AIMAG(UL(J1))) 

J1=J1+N 

IF(RN-ULM) 1,2,2 

1 RN=ULM 

2 CONTINUE 
SCL(1)=1./RN 

5 CONTINUE 
NM1=N-1 
K2=0 

DO 17 K=1,NM1 
BIG=0. 

DO 11 l=K,N 
IP=IPS(I) 

IPK=IP+K2 

SIZE=(ABS(REAL(UL(IPK)))+ABS(AIMAG(UL(IPK))))*SCL( 

IP) 

IF(SIZE-BIG) 11,11,10 

10 BIG=SIZE 
IPV=I 

11 CONTINUE 
IF(IPV-K) 14,15,14 

14 J=IPS(K) 

IPS(K)=IPS(IPV) 

IPS(IPV)=J 

15 KPP=IPS(K)+K2 
PIVOT=UL(KPP) 

KP1=K+1 

DO 16 l=KP1,N 
KP=KPP 
IP=IPS(I)+K2 
EM=-UL(IP)/PlVOT 
18 UL(IP)=-EM 
DO 16 J=KP1,N 
IP=IP+N 
KP=KP+N 

UL(IP)=UL(IP)+EM*UL(KP) 

16 CONTINUE 
K2=K2+N 

17 CONTINUE 



RETURN 

END 

FUNCTION BLOG(X) 

IF(X.GT.,1) GOTO 1 
X2=X*X 

BLOG=((.075*X2-.1666667)*X2+1 .)*X 
RETURN 

1 BLOG=ALOG(X+SQRT(1 ,+X*X)) 

RETURN 

END 

SUBROUTINE ZLOAD(NP,RH,ZHZ0,Z) 

C 

C COMPUTES IMPEDANCE MATRIX ELEMENTS FOR 
LOADED BODIES OF REV 
C Z0(l) IS THE SURF IMPEDANCE OF THE ITH 
SEGMENT (NP-1 SEGMENTS) 

C Z(.) ARE THE IMPEDANCE MATRIX TERMS 
(TRIDIAGONAL FOR T-T 

C SUBMATRIX; DIAGONAL FOR P-P SUBMATRIX). 
STORED IN COL VECTOR. 

C 

COMPLEX 

Cl ,C2,Z0(400),Z(1500),X1 ,X2,X3.V1 ,Y2,Y3,FN(400) 
COMPLEX U1.U2.U3, XI, Yl 

DIMENSION RH(400).ZH(400),RS(400),D(400),SV(400) 

Pl=3. 14159 

MT=NP-2 

MP=NP*1 

N=MT+MP 

DO 10 IP=2,NP 

ll=IP-1 

DR=RH(IP)-RH(II) 

DZ=ZH(IP)-ZH(II) 

D(II)=SQRT(DR*DR+DZ‘DZ) 

SV(II)=DR/D(II) 

RS(II)=.5*(RH(IP)+RH(II)) 

DS=D(ll)*SV(ll)/2. 

Q1=RS(II)+DS 
Q2=RS(1I)-DS 
FN(II)=1 . 

IF((ABS(Q2).GT.1 .E*6).AND.(ABS(Q1).GT.1 .E*6)) 

* FN(ll)=ALOG(Q1/Q2) 

10 CONTINUE 
L0=MT*3-2 
DO 20 1=1, MP 
C1=PI*Z0(I) 

IF(I.EQ.MP) GOTO 80 
Kl=2 

IF(I.EQ.1 ) Kl=1 
IF(I.EQ.MT) Kl=3 
11 = 1+1 

C2=PI*Z0(II) 

A=SV(I) 

IF(ABS(A).LT.1 E-6) GO TO 41 

X1=C1*FN(l)/2./A 

X2=C1 *2./A*(1 .*RS(I)*FN(I)/D(I)/A) 

X3=*X2*RS(I)/D(I)/A 
GO TO 42 

41 CONTINUE 
X1=C1/2./RS(l)*D(l) 

X2=(0.,0.) 

X3=C1*D(l)/6./RS(l) 

42 CONTINUE 
A=SV(ll) 

IF(ABS(A).LT.1 .E*6) GO TO 45 
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Y1=C2TN(ll)/2./A 

Y2=C2*2./A*(1 .-RS(II)*FN(II)/D(II)/A) 
Y3=-Y2*RS(II)/D(II)/A 
GO TO 40 
45 CONTINUE 
Yl=C2/2./RS(ll)*D(ll) 

Y2=(0.,0.) 

Y3=C2*D(ll)/6./RS(ll) 

40 CONTINUE 
C 

C DEFINE TRIDIAGONAL ELEMENTS FOR T*T 
SUBMATRIX (STORED IN COLS) 

C (U1- DIAG; U2- LOWER; U3- UPPER) 

C 

XI=X1+X2+X3 
YI=Y1-Y2+Y3 
IF(KI.EQ.I) XI=C1/SV(I) 

C IF(KI.EQ.3) YI=C2/SV(II) 

U1=XI+YI 

U2=X1-X3 

U3=Y1-Y3 

L=2+(l-2)*3 

IF(KI.EQ.I) L=0 

L1=L+1 

L2=L+2 

L3=L+3 

goto (50, 60,70), ki 
50 Z(L1)=U1 
Z(L2)=U2 
GOTO 80 
60 Z(L1)=U3 
Z(L2)=U1 
Z(L3)=U2 
GO TO 80 
70 Z(L1)=U3 
Z(L2)=U1 

80 Z(L0+I)=2.*C1 *D(I)/RS(I) 

20 CONTINUE 
RETURN 
END 

SUBROUTINE ZTOT(MT,MP,ZL.Z) 

C 

C ADDS THE SURF IMPEDANCE TERMS TO THE 
TRIDIAGONAL ELEMENTS OF 
C THE BOR IMPEDANCE MATRIX Z. 

C 

COMPLEX ZL(700),Z(50000) 

N=MT+MP 
M0=MT*3-2 
DO 100 1=1, MP 
L0=MT*N+(M)*N+MT 
IF(I.EQ.MP) GOTO 80 
Kl=2 

IF(I.EQ.I) Kl=1 

IF(I.EQ.MT) Kl=3 

L2=(I-1)*N+I 

L1=L2-1 

L3=L2+1 

M=2+3*(l*2) 

IF(KI.EQ.I) M=0 
M1=M+1 
M2=M+2 
M3=M+3 

goto (50, 60,70), ki 
50 Z(L2)=Z(L2)+ZL(M1) 



Z(L3)=Z(L3)+ZL(M2) 

GO TO 80 

60 Z(L1)=Z(L1)+ZL(M1) 

Z(L2)=Z(L2)+ZL(M2) 

Z(L3)=Z(L3)+ZL(M3) 

GO TO 80 

70 Z(L1)=Z(L1)+ZL(M1) 

Z(L2)=Z(L2)+ZL(M2) 

80 Z(L0+l)=Z(L0+l)+ZL(M0+l) 

100 CONTINUE 
RETURN 
END 

subroutine interp(u,v,x,y,nn) 

c program to interpolate linearly between arrays of x and y 
c values (with dimensions nn). u is specified and v is 
returned 

c uniform sampling in x (the ordinate) is assumed, 
dimension x(nn),y(nn) 
delx=x(2)-x(1) 
n1=int(u/delx)+1 
n2=n1+1 
dely=y(n2)-y(n1) 
sgn=sign(1 ..dely) 
dely=abs(dely) 
alpha=u-x(n1) 

2 eta=dely*alpha'detx 

v=y(n1)+sgn*zeta 

return 

end 

SUBROUTINE 

PLANEW(NWIRES,NW1,NW2,NP1,NP2,XH,YH,ZH. 

* THR,PHR,R) 

C 

C PLANE WAVE EXCITATION VECTOR ELEMENTS FOR 
WIRE AND 

C ATTACHMENT SEGMENT. INCIDENCE DIRECTION IS 
(THR.PHR). 

C - ONLY ONE INCIDENCE ANGLE PER CALL 

C - WIRE DOES NOT HAVE TO LIE IN THE Z=0 PLANE 

C 

COMPLEX 

U0,AA,BB,C,R(400),CEXP,EXP.FI1,FI2,SI,DI,CMPLX 

DIMENSION 

ZH(400),NS(4),NW1(4),NW2(4),XH(400),YH(400) 
MP2=NP2*1 
MT2=NP2-2 
MT22=2*MT2 
DO 5 L=1 ,NWIRES 
5 NS(L)=NW2(L)*NW1(L)+1 
U0=(0.,0.) 

CC=COS(THR) 

SS=SIN(THR) 

CP=COS(PHR) 

SP=SIN(PHR) 

UP=SS*CP 

VP=SS*SP 

WP=CC 

DO 12 IP=1 ,MP2 

ll=IP+NP1 

1 = 11+1 

ZS=.5*(ZH(I)+ZH(II)) 

XS=.5*(XH(I)+XH(II)) 

YS=.5*(YH(I)+YH(II)) 

DX=XH(I)-XH(II) 

DY=YH(I)-YH(II) 



58 



ooo ooo ooo ooo 



D1=SQRT(DX**2+DY**2) 

SU=DY/D1 

CU=DX/D1 

C FOR WIRES IN PERPENDICULAR TO 2 SIN(V)=1 AND 
COS(V)=0 
SV=1.0 
CV=0.0 

WIRE SEGMENT CALCULATIONS 

A=UP*CU+VP*SU 

B=UP*XS+VP*YS+WP*ZS 

C=CMPLX(0.,A) 

EXP=CEXP(CMPLX(0.,B)) 

AA=CC*(CU*CP+SU*SP) 

BB=SU*CP-SP*CU 

PSI=D1*A/2. 

IF(PSI.NE.O.) GO TO 60 
SINC=1. 

GO TO 61 

60 SINC=SIN(PSI)/PSI 

61 COSP=COS(PSI) 

FI1=SINC*D1*EXP/2. 

FI2=(0.,0.) 

IF(ABS(A).LT.1 .E-4) GO TO 62 
CSP=COSP-SINC 
IF(ABS(CSP).LT.1.E-4) GO TO 62 
FI2=EXP/C*CSP 

62 CONTINUE 
SI=FI1+FI2 
DI=F11-FI2 

R-WIRE-THETA 

IF(IP.EQ.MP2) GOTO 10 
R(IP)=AA*SI 
R(IP+MT2)=BB*SI 
10 CONTINUE 

R-WIRE-PHI 

14 IF(IP.EQ.I) GO TO 12 
R(IP-1)=R(IP-1)+AA*DI 
R(1P-1+MT2)=R(IP-1+MT2)+BB*DI 
12 CONTINUE 

WRITE OVER DUMMY WIRE SEGMENTS 

c IF(NWIRES.EQ.I) GO TO 210 
c NSPTS=0 
c DO 209 L=1 ,NWIRES-1 
c NSPTS=NSPTS+NS(L) 
c IN=NSPTS-2 
c R(IN+1 )=U0 
c R(IN+2)=U0 
c IN=IN+MT2 
c R(IN+1)=U0 
c R(IN+2)=U0 
c 209 CONTINUE 
210 RETURN 
END 

SUBROUTINE 

PLANES(M1 > NP0 t NP1,NP2,NT,RH.ZH,XT,AT,THR,R) 

C 



C PLANE WAVE EXCITATION VECTOR FOR BOR 
ELEMENTS. 

C NO CHANGE FROM HARRINGTON'S ORIGINAL 
PROGRAM OTHER THAN 
C 

C ** ONLY ONE MODE PER CALL ** 

C 

COMPLEX 

R(1000),U,U1 ,UA 1 UB,FA(50),FB(50) 1 F2A,F2B,F1 A,F1 B 
COMPLEX U2,U3,U4,U5,CMPLX 
DIMENSION RH(400),ZH(400),XT(4),AT(4) 
DIMENSION R2(4) r 22(4),BJ(5000) 

U=(0. T 1.) 

U 1=3.1 41 593*U**M1 

NF=1 

M2=M1 

NP=NP0+NP1+NP2 

NRF=NP0+NP1 

MP=NP-1 

MT=MP-1 

N=MT+MP 

N2=2*N 

CC=COS(THR) 

SS=SIN(THR) 

M3=M1 + 1 

M4=M2+3 

IF(MI.EQ.O) M3=2 

M5=M1+2 

M6=M2+2 

DO 12 IP=1 ,MP 

K2=IP 

l=IP+1 

DR=.5*(RH(I)-RH(IP)) 

DZ=.5*(ZH(I)-ZH(IP)) 

D1=SQRT(DR*DR+DZ‘DZ) 

R1=.25*(RH(I)+RH(IP)) 

IF(ABS(R1).LT.1.E-5) R1 = 1. 

Z1=.5*(ZH(I)+ZH(IP)) 

DR=.5*DR 
D2=DR/R1 
DO 13 L=1 ,NT 
R2(L)=R1+DR*XT(L) 

Z2(L)=Z1+DZ*XT(L) 

13 CONTINUE 
D3=DR*CC 
D4=-DZ*SS 
D5=D1*CC 
DO 23 M=M3,M4 
FA(M)=0. 

FB(M)=0. 

23 CONTINUE 
DO 15 L=1 ,NT 
X=SS*R2(L) 

IF(X.GT..5E-7) GOTO 19 
DO 20 M=M3,M4 
BJ(M)=0. 

20 CONTINUE 
BJ(2)=1 . 

S=1 . 

GOTO 18 
19 M=2.8*X+14.-2./X 
IFPC.LT..5) M=1 1 .8+ALOG 1 0(X) 

IF(M.LT.M4) M=M4 
BJ(M)=0. 

JM=M-1 
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BJ(JM)=1 . 

DO 16 J=4,M 
J2=JM 
JM=JM-1 
J1=JM-1 

BJ(JM)=J1/X*BJ(J2)-BJ(JM+2) 

16 CONTINUE 
S=0. 

IF(M.LE.4) GOTO 24 
DO 17 J=4,M,2 
S=S+BJ(J) 

17 CONTINUE 

24 S=BJ(2)+2.*S 

18 ARG=Z2(L)*CC 

UA=AT(L)/S*CMPLX(COS(ARG),SIN(ARG)) 
UB=XT(L)*UA 
DO 25 M=M3,M4 
FA(M)=BJ(M)*UA+FA(M) 

FB(M)=BJ(M)’*UB+FB(M) 

25 CONTINUE 
15 CONTINUE 

IF(MI.NE.O) GO TO 26 
FA(1)=- FA(3) 

FB(1)=- FB(3) 

26 UA=U1 

DO 27 M=M5.M6 

M7=M-1 

M8=M+1 

F2A=UA’*(FA(M8)+FA(M7)) 

F2B=UA*(FB(M8)+FB(M7)) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7)) 

F1B=UB*(FB(M8)-FB(M7)) 

U4=D4*UA 

U2=D3*F1 A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT)=-D5*(F2A+D2*F2B) 

R(K5+MT)=D1 *(F1A+D2*F1B) 

IF(IP.EQ.I) GO TO 21 
R(K1)=R(K1)+U2-U3 
R(K4)=R(K4)+U4-U5 
IF(IP.EQ.MP) GO TO 22 

21 R(K2)=U2+U3 
R(K5)=U4+U5 

22 K2=K2+N2 
UA=UB 

27 CONTINUE 
12 CONTINUE 

c R(NP0-1)=(0.,0.) 
c R(NPO)=(0.,0.) 
c R(NP0+MT)=(0.,0.) 
c R(NRF-1)=(0.,0.) 
c R(NRF)=(0.,0.) 
c R(NRF+MT)=(0.,0.) 

99 RETURN 
END 

SUBROUTINE 

ZMATSS(M1 ,M2,NP0,NP1 ,NP2,NPHI,NT,RH,ZH,X,A,XT,A 
Z) 



C THREE DETACHED SURFACES 

C SURFACE-SURFACE IMPEDANCE ELEMENTS. 
REMAINS UNCHANGED FROM HARRINGTON 
C EXCEPT MULTIPLE SURFACES ARE PERMITTED. 

THE FIRST SURFACE HAS NPO 

C POINTS, THE SECOND NP1 POINTS AND THE THIRD 
NP2. 

C NPO PARABOLOID PTS; NP1 HYPERBOLOID PTS; 

NP2 CAVITY PTS. 

C 

COMPLEX 

Z(50000),U1 ,U2,U3,U4,U5,U6,U7,U8,U9,UA,UB,G4A(4),G5 
A(4) 

COMPLEX 

G6A(4),G4B(4),G5B(4),G6B(4),H4A,H5A,H6A,H4B,H5B 
COMPLEX CMPLX,H6B,UC,UD.GA(400),GB(400) 
DIMENSION 

RH(400)ZH(400),X(400),A(400),AT(4),RS(400),ZS(400) 

DIMENSION 

D(400),DR(400),DZ(400),DM(400),C2(400),C3(400),R2(4) 

DIMENSION 

C4(400),C5(400),C6(400),Z7(4),R7(4),Z8(4),R8(4) 
DIMENSION XT(4),Z2(4) 

CT=2. 

CP=.1 

NP=NP0+NP1+NP2 
NRF=NP0+NP1 
DO 10 1=2, NP 
12 = 1-1 

RS(I2)=.5*(RH(I)+RH(I2)) 

ZS(I2)=.5*(ZH(I)+ZH(I2)) 

D1=.5*(RH(I)-RH(I2)) 

D2=.5*(ZH(I)-ZH(I2)) 

D(I2)=SQRT(D1*D1+D2*D2) 

DR(I2)=D1 

DZ(I2)=D2 

DM(I2)=D(I2)/RS(I2) 

10 CONTINUE 
M3=M2-M1+1 
M4=M1-1 
PI2=1. 570796 
DO 11 K=1 ,NPHI 
PH=PI2’(X(K)+1 .) 

C2(K)=PH*PH 

SN=SIN(.5*PH) 

C3(K)=4.*SN*SN 

A1=PI2*A(K) 

D4=.5*A1 *C3(K) 

D5=A1*COS(PH) 

D6=A1*SIN(PH) 

M5=K 

DO 29 M=1 ,M3 

PHM=(M4+M)*PH 

A2=COS(PHM) 

C4(M5)=D4*A2 

C5(M5)=D5*A2 

C6(M5)=D6*SIN(PHM) 

M5=M5+NPHI 
29 CONTINUE 

11 CONTINUE 
MP=NP-1 
MT=MP-1 
N=MT+MP 
N2N=MT’N 
N2=N*N 
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U1=(0.,.5) 

U2=(0.,2.) 

JN=-1-N 

DO 15 JQ=1,MP 

KQ=2 

IF(JQ.EQ.I) KQ=1 
IF(JQ.EQ.MP) KQ=3 
R1=RS(JQ) 

21=ZS(JQ) 

D1=D(JQ) 

D2=DR(JQ) 

D3=DZ(JQ) 

D4=D2/R1 

D5=DM(JQ) 

SV=D2/D1 

CV=D3/D1 

T6=CT*D1 

T62=T6+D1 

T62=T62*T62 

R6=CP*R1 

R62=R6*R6 

DO 12 L=1 ,NT 

R2(L)-R1+D2*XT(L) 

Z2(L)=Z1+D3*XT(L) 

12 CONTINUE 
U3=D2*U1 
U4=D3*U1 
DO 16 IP=1 ,MP 
R3=RS(IP) 

Z3=ZS(IP) 

R4=R1-R3 

Z4=Z1-Z3 

FM=R4*SV+Z4*CV 

PHM=ABS(FM) 

PH=ABS(R4*CV-Z4*SV) 

D6=PH 

IF(PHM.LE.D1 ) GO TO 26 
D6=PHM-D1 

D6=SQRT(D6*D6+PH*PH) 

26 IF(IP.EQ.JQ) GO TO 27 
KP=1 

IF(T6.GT.D6) KP=2 
IF(R6.GT.D6) KP=3 
GO TO 28 

27 KP=4 

28 GOTO (41 ,42,41 ,42),KP 
42 DO 40 L=1 ,NT 

D7-R2(L)*R3 

D8=Z2(L)-Z3 

Z7(L)=D7*D7+D8*D8 

R7(L)=R3*R2(L) 

Z8(L)=.25*Z7(L) 

R8(L)=.25*R7(L) 

40 CONTINUE 
Z4=R4*R4+Z4*Z4 
R4=R3*R1 
R5=.5*R3*SV 
DO 33 K=1 ,NPHI 
A1=C3(K) 

RR=Z4+R4*A1 

UA=0. 

UB=0. 

IF(RR.LT.T62) GO TO 34 
DO 35 L=1,NT 
R=SQRT(Z7(L)+R7(L)*A1) 



SN=*SIN(R) 

CS=COS(R) 

UC=AT(L)/R*CMPLX(CS,SN) 

UA=UA+UC 

UB=XT(L)*UC+UB 

35 CONTINUE 
GO TO 36 

34 DO 37 L=1,NT 
R=SQRT(Z8(L)+R8(L)*A1 ) 
SN=-SIN(R) 

CS=COS(R) 

UC=AT(L)/R*SN*CMPLX(*SN,CS) 

UA=UA+UC 

UB=XT(L)*UC+UB 

37 CONTINUE 
A2=FM+R5*A1 
D9=RR-A2*A2 
R=ABS(A2) 

D7=R-D1 

D8=R+D1 

D6=SQRT(D8*D8+D9) 

R=SQRT (D7*D7+D9) 

IF(D7.GE.O.) GO TO 38 
A1=ALOG((D8+D6)’(-D7+R)/D9) / D1 
GO TO 39 

38 A1=ALOG((D8+D6)/(D7+R))/D1 

39 UA=A1+UA 

UB=A2*(4./(D6+R)‘A1)/D1+UB 

36 GA(K)=UA 
GB(K)=UB 

33 CONTINUE 
K1=0 

DO 45 M=1 ,M3 
H4A=0. 

H5A=0. 

H6A=0. 

H4B=0. 

H5B=0. 

H6B=0. 

DO 46 K=1 ,NPHI 

K1=K1+1 

D6=C4(K1) 

D7=C5(K1) 

D8=C6(K1 ) 

UA=GA(K) 

UB=GB(K) 

H4A=D6*UA+H4A 
H5A=D7*UA+H5A 
H6A=D8*UA+H6A 
H4B=D6*UB+H4B 
H5B=D7*UB+H5B 
H6B=D8*UB+H6B 
46 CONTINUE 
G4A(M)=H4A 
G5A(M)=H5A 
G6A(M)=H6A 
G4B(M)=H4B 
G5B(M)=H5B 
G6B(M)=H6B 
45 CONTINUE 
IF(KP.NE.4) GO TO 47 
A2=D1/(PI2*R1) 

D6=2./D1 

D8=0. 

DO 63 K=1 ,NPHI 
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A1=R4*C2(K) 

R=R4*C3(K) 

IF(R.LT.T62) GO TO 64 
D7=0. 

DO 65 L=1 ,NT 

D7=D7+AT(L)/SQRT(27(L)+A1) 

65 CONTINUE 
GO TO 66 

64 A1=A2/(X(K)+1 .) 

D7=D6’ALOG(AUSQRT(1 .+A1 *A1)) 

66 D8=D8+A(K)*D7 
63 CONTINUE 

A1=.5*A2 

A2=1./A1 

D8=-PI2’D8+2./R 1 *(BLOG(A2)+A2*BLOG(A1 )) 

DO 67 M=1 ,M3 
G5A(M)=D8+G5A(M) 

67 CONTINUE 
GO TO 47 

41 DO 25 M=1.M3 
G4A(M)=0. 

G5A(M)=0. 

G6A(M)=0. 

G4B(M)=0. 

G5B(M)=0. 

G6B(M)=0. 

25 CONTINUE 
DO 13 L=1,NT 
A1=R2(L) 

R4=A1-R3 

Z4=Z2(L)-Z3 

Z4=R4*R4+Z4*Z4 

R4=R3*A1 

DO 17 K=1 ,NPHI 

R=SGRT(Z4+R4*C3(K)) 

SN=-SIN(R) 

CS=COS(R) 

GA(K)=CMPLX(CS,SN)/R 
17 CONTINUE 
D6=0. 

IF(R62.LE.Z4) GO TO 51 
DO 62 K=1 ,NPHI 

D6=D6+A(K)/SQRT(Z4+R4*C2(K)) 

62 CONTINUE 
Z4=3 . 1 4 1 593/SQRT (Z4/R4) 

D6=-PI2*D6+ALOG(Z4+SQRT(1 .+Z4’Z4)) SQRT(R4) 
51 A1=AT(L) 

A2=XT(L)*A1 

K1=0 

DO 30 M=1 ,M3 
U5=0. 

U6=0. 

U7=0. 

DO 32 K=1,NPHI 
UA=GA(K) 

K1=KH1 

U5=C4(K1)’UA+U5 
U6=C5(K1)*UA+U6 
U7=C6(K1)*UA+U7 
32 CONTINUE 
U6ssD6+U6 

G4A(M)=A1 *U5+G4A(M) 

G5A{M)=AVU6+G5A(M) 

G6A(M)=A1*U7+G6A(M) 

G4B(M)=A2*U5+G4B(M) 



G5B(M)=A2*U6+G5B(M) 

G6B(M)=A2*U7+G6B(M) 

30 CONTINUE 

13 CONTINUE 

47 A1=DR(IP) 

UA=A1*U3 

UB=DZ(IP)*U4 

A2=D(IP) 

D6=-A2*D2 

D7=D1*A1 

D8=D1*A2 

JM=JN 

DO 31 M=1,M3 

FM=M4+M 

A1=FM*DM(IP) 

H5A=G5A(M) 

H5B=G5B(M) 

H4A=G4A(M)+H5A 

H4B=G4B(M)+H5B 

H6A=G6A(M) 

H6B=G6B(M) 

U7=UA*H5A+UB*H4A 

U8=UA*H5B+UB*H4B 

U5=U7*U8 

U6=U7+U8 

U7=-U1*H4A 

U8=D6*H6A 

U9=D6*H6B-A1 *H4A 

UC=D7 4 (H6A+D4*H6B) 

UD=FM*D5*H4A 

K1=IP+JM 

K2=K1+1 

K3=K1+N 

K4=K2+N 

K5=K2+MT 

K6=K4+MT 

K7=K3+N2N 

K8=K4+N2N 

GO TO (18, 20, 19), KQ 

18 Z(K6)=U8+U9 
IF(IP.EQ.I) GO TO 21 
Z(K3)=Z(K3)+U6*U7 
Z(K7)=Z(K7)+UC*UD 
IF(IP.EQ.MP) GO TO 22 

21 Z(K4)=U6+U7 
Z(K8)=UC+UD 
GO TO 22 

19 Z(K5)=Z(K5)+U8-U9 
IF(IP.EQ.I) GO TO 23 
Z(K1 )=Z(K1 )+U5+U7 
Z(K7)=Z(K7)+UC*UD 
IF(IP.EQ.MP) GO TO 22 

23 Z(K2)=Z(K2)+U5-U7 
Z(K8)=UC+UD 

GO TO 22 

20 Z(K5)=Z(K5)+U8-U9 
Z(K6)=U8+U9 
IF(IP.EQ.I) GOTO 24 
Z(K1 )=Z(K1 )+U5+U7 
Z(K3)=Z(K3)+U6-U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 

24 Z(K2)=Z(K2)+U5-U7 
Z(K4)=U6+U7 
Z(K8)=UC+UD 
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22 2(K8+MT)=U2*(D8*(H5A+D4*H5B)-A1*UD) 
JMssJM+N2 
31 CONTINUE 
16 CONTINUE 
JN=JN+N 
15 CONTINUE 
C 

C THREE MULTIPLE BODIES USING THE SIMPLIFIED 
APPROACH 

C NULL OUT ROWS AND COLS FOR THE FIRST 
DUMMY SEGMENT 
C 

DO 100 LSS=1 ,N 
LS=LSS-1 

Z(LS*N+NP0)=(0.,0.) 

Z(LS*N+NP0-1)=(0.,0.) 

100 Z(LS*N+NP0+MT)=(0.,0.) 

DO 101 LS^I.N 
Z((NP0*2)*N+LS)=(0.,0.) 

Z((NP0-1)*N+LS)=(0.,0.) 

101 Z((NP0-1+MT)*N+LS)=(0.,0.) 
Z((NP0-2)*N+NP0-1)=(1 .,0.) 

Z((NP0-1 )*N+NP0)=(1.,0.) 
Z((MT+NP0-1)*N+NP0+MT)=(1.,0.) 

C 

C NULL OUT ROWS AND COLS FOR THE SECOND 
DUMMY SEGMENT 
C 

DO 200 LSS=1 ,N 
LS=LSS-1 

Z(LS*N+NRF)=(0.,0.) 

Z(LS*N+NRF-1)=(0.,0.) 

200 Z(LS*N+NRF+MT)=(0.,0.) 

DO 201 LS=1 ,N 
Z((NRF-2)*N+LS)=(0.,0.) 

Z((NRF-1)*N+LS)=(0.,0.) 

201 Z((NRF-1+MT)*N+LS)=(0.,0.) 
Z((NRF-2)*N+NRF-1)=(1.,0.) 
Z((NRF-1)*N+NRF)=(1.,0.) 
Z((MT+NRF-1)*N+NRF+MT)=(1.,0.) 

999 RETURN 
END 

SUBROUTINE ZASMB0(NP1 ,NP2. MODES, Z,ZWW) 
C 

C ASSEMBLES MODE INDEPENDENT BLOCKS OF Z 
C 

COMPLEX Z(500000),ZWW( 10000) 

MT1=NPl-2 
MP1=NP1-1 
N=MT1+MP1 
MT2=NP2-2 
NBLOCK=2*MODES+1 
NS=NBLOCK*N 
NR=NS+MT2 
NS2=NS*NR 
NWW=NS2+NS 
DO 10 IC=1 ,MT2 
IC1=IC-1 
DO 10 IR=1 ,MT2 
K=NWW+IC1*NR+IR 
IW=IC1*MT2+IR 
10 Z(K)*ZWW(IW) 

RETURN 

END 



SUBROUTINE 

ZASMBN(NP1,NP2,MODES,NM,Z.ZSS2SW) 

C 

C ASSEMBLES MODE DEPENDENT BLOCKS IN Z. 
MODES=TL # OF MODES; 

C NM IS THE MODE # OF THE BLOCKS ZSS.ZSW.ZSJ. 
(IF NM<0 THEN 

C ZSS IS OMITTED - ONLY FILLED FOR POSITIVE NM) 
C 

COMPLEX ZSS(50000),Z(500000),ZSW(50000) 

MT1=NPl-2 

MPlsNPM 

MT2=NP2-2 

N=MP1+MT 1 

NBLOCK=2*MODES+1 

NS=NBLOCK*N 

NR=NS+MT2 

MT=MT1+MT2 

MTN=MT*N 

IF(NM.LT.O) GOTO 225 
C 

C ZSS BLOCKS FOR MODE NM (+ AND * FILLED IN 1 
LOOP) 

C 

NTOP1=MODES-NM 
NTOP2=NBLOCK-(NTOP1 + 1 ) 
NSS2=NTOP1*(N+NR*N) 

NSS1=NTOP2*(N+NR*N) 

DO 200 IC=1 ,MP1 
IC1=IC-1 

DO 200 IR=1,MP1 
IF(IC.GT.MTI) GOTO 100 
IF(IR.GT.MTI) GOTO 50 
I0=IC1*N+IR 
ITT=IC1*NR+IR 
Z(ITT+NSS1)=ZSS(I0) 

IF(NM.NE.O) Z(ITT+NSS2)=ZSS(I0) 

50 I0=MT1+ICTN+IR 
IPT=ICTNR+MT1 + IR 
Z(IPT+NSS1)= ZSS(IO) 

IF(NM.NE.O) Z(IPT+NSS2)=-ZSS(I0) 

100 IF(IR.GT.MTI) GOTO 150 
I0=N*MT 1+ICTN+IR 
ITP=MT 1 ‘NR+IC1 *NR+1R 
Z(ITP+NSS1)= ZSS(IO) 

IF(NM.NE.O) Z(ITP+NSS2)=-ZSS(I0) 

150 I0=N*MT1+MT1 + ICTN+IR 
IPP=1C1 *NR+MT1 +MT1 *NR+IR 
Z(IPP+NSS1)=ZSS(I0) 

IF(NM.NE.O) Z(IPP+NSS2)=ZSS(I0) 

200 CONTINUE 
225 CONTINUE 
C 

C ZSW AND ZWS BLOCKS 
C 

NTOP1=MODES+NM 

NTOP2=NBLOCK- (NTOP1 + 1 ) 

NSW=NTOPTN+NS*NR 

NWS=NTOP2*NR*N+NS 

DO 250 IC=1 ,MT2 

IC1=IC-1 

DO 250 IR=1 ,N 

ISW=IR+IC1*NR+NSW 

IWS=(IR-1)*NR+IC+NWS 

I0=ICTN+IR 
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Z(ISW)=ZSW(IO) 

Z(IWS)=ZSW(IO) 

250 CONTINUE 
RETURN 
END 

SUBROUTINE 

ZMATSW(NWIRES,NW1,NW2,NP11,NP12,NP13 I 

•XH,YH,RH,ZH,NT,XT,AT,NG,XG,AG.MODE,AZ) 

C 

c »»»»»> revised version ««««« 
c 

C CALCULATION OF THE SURFACE-WIRE MATRIX 
ELEMENTS. ATTACHMENT 

C POINT NOT AT A SURFACE EDGE, BUT CAN BE ON Z 
AXIS (IC#1). SURFACE 

C START INDEX IS 1. TWO BORS OF LENGTHS NP11, 
NP12AND NP13. WIRE 

C START AND STOP POINTS ARE NW1 AND NW2. 

C 

COMPLEX 

CEXP,EXP,Z(50000),G5,G6,G7,CONH.F5,F6,F7 

COMPLEX 

U,U0,PSI,T0,T 1 ,T2,T3,ST1 ,ST2,ST3,CON,SPl ,SP2,SP3 
COMPLEX U5,U6,U7,U8 t U9,CMPLX,T4.T5,T6,T7 
DIMENSION 

RS1(400),RH(400),ZH(400),XG(400),AG(400),XT(4),AT(4) 

DIMENSION 

ZS1 (400),D1 (400),D(400),S(400),C(400),RS(400) 
DIMENSION 

ZS(400),S1(400),C1(400),NW1(4),NW2(4).NS(4) 

DIMENSION 

XH(400),YH(400),XS1 (400),YS1 (400),CU(400),SU(400) 
Pl=3, 14159 
PI2=2.*PI 

NP1=NP1 1+NP12+NP13 
NRF=NP1 1+NP12 
DO 1 L=1 .NWIRES 
1 NS(L)=NW2(L)-NW1(L)+1 
U0=(0.,0.) 

U=(0.,1.) 

c 

C DEFINE SURFACE GEOMETRY TERMS 
C 

MP=NP1-1 
MT=MP-1 
NROW=MP+MT 
DO 5 1=2, NP1 
12 = 1-1 

RS(I2)=.5*(RH(I)+RH(I2)) 

ZS(I2)=.5’(ZH(I)+ZH(I2)) 

DR=RH(I)-RH(I2) 

DZ=ZH(I)-ZH(I2) 

D(I2)=SQRT(DR**2+DZ**2) 

S(I2)=DR/D(I2) 

5 C(I2)=DZ/D(I2) 

C 

C DEFINE WIRE GEOMETRY TERMS 
C 

NW=NW2(NWIRES)-NW1 (1 ) 

NWP=NW+1 

MT2=NWP-2 

DO 10 N=2,NWP 

N0=N-1 

|=NP1+N 

12 = 1-1 



RS1(N0)=.5*(RH(I)+RH(I2)) 

ZS1(N0)=.5*(ZH(I)+ZH(I2)) 

XS1(N0)=.5*(XH(I)+XH(I2)) 

YS1 (N0)=.5*(YH(I)+YH(I2)) 

DX=XH(I)-XH(I2) 

DY=YH(I)-YH(I2) 

D1(N0)=SQRT(DX**2+DY**2) 

C FOR WIRES PARALLEL TO THE XY PLANE SIN(V)=1 
AND COS(V)=0 
S1(N0)=1. 

Cl (N0)=0. 

CU(N0)=DX/D1 (NO) 

SU(N0)=DY/D1 (NO) 

10 CONTINUE 
CON=U/PI2/2. 

C INTEGRATING FROM 0 TO 2*PI 
01=PI 
02=PI 

JN=-1-NROW 

C 

C SEGMENT LOOPS: q=WIRE, p=SURFACE 
C 

DO 100 JQ=1,NW 
KQ=2 

IF(JQ.EQ.I) KQ=1 
IF(JQ.EQ.NW) KQ=3 
H1=D1(JQ)/2. 

ZW=ZS1(JQ) 

DO 22 IP=1 ,MP 

ST1=U0 

ST2=U0 

ST3=U0 

SP1=U0 

SP2=U0 

SP3=U0 

AA=S(IP)*CU(JQ) 

BB=S(IP)*SU(JQ) 

RD=RS(IP)*D1(JQ) 

DD=D(IP)*D1(JQ) 

CONH=CON*D(IP) 

C 

C FIRST TERM (H- INTEGRATION) 

C 

DO 30 1=1, NT 
H=H1*XT(I) 

HHD=2.’H D1(JQ) 

XW=XS1(JQ)+H*CU(JQ) 

YW=YS1 (JQ)+H’SU(JQ) 

C 

C PHI INTEGRATIONS FOR G-HAT FUNCTIONS 
C 

G7=U0 

G5=U0 

G6=U0 

DO 32 K=1 ,NG 

PH=01*XG(K)+02 

PHM=PH'FLOAT(MODE) 

XS=RS(lP)*COS(PH) 

YS=RS(IP)*SIN(PH) 

EXP=CEXP(CMPLX(0.,-PHM)) 

RP=SQRT((XS-XW)”2+(YS-YW)’*2+(ZS(IP)-ZW)**2) 

PSI=CEXP(CMPLX(0.,-RP))/RP 

F5=EXP*COS(PH)*PSI 

F6=EXP*SIN(PH)*PSI 

F7=EXP’PSI 
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G7=G7+F7*AG(K) 

G5=G5+F5*AG(K) 

32 G6=G6+F6*AG(K) 

G7-01*G7 

G5=01*G5 

G6=01'G6 

T0=G5*AA/4. 

T1=T0*HHD 

T2=G6*BB/4. 

T3=T2*HHD 
T 4=G5*SU(JQ)/2. 

T5=T4*HHD 

T6=*G6*CU(JQ)/2. 

T7=T6*HHD 

ST1=ST1+AT(I)*(T0+T2) 

ST2=ST2+AT(I)*(T 1 +T3) 

ST3=ST3+AT(I)'G7/DD 

SP1=SP1+AT(I)*(T4+T6) 

SP2=SP2+AT(I)*(T5+T7) 

SP3=SP3+AT(l)*U*G7*FLOAT(MODE)/RD 

30 CONTINUE 
C 

C ZSW-t TERMS (K1 ,K2,K3,K4) 

C 

U5=(ST1-ST2)*CONH*H1 

U6=(STl+ST2)*CONH*H1 

U7=-ST3*CONH*Hl 

C 

C ZSW-p TERMS (K5,K6) 

C 

U8=SP1 *CONH*H1 

U9=(SP2+SP3)*CONH*Hl 

K1=IP+JN 

K2=K1+1 

K3=K1+NROW 

K4=K2+NROW 

K5=K2+MT 

K6=K4+MT 

GO TO (18, 20,19), KQ 

1 8 Z(K6)=U8+U9 
IF(IP.EQ.I) GO TO 21 
Z(K3)=Z(K3)+U6-U7 
IF(IP.EQ.MP) GO TO 22 

21 Z(K4)=U6+U7 
GO TO 22 

19 Z(K5)=Z(K5)+U8-U9 
IF(IP.EQ.I) GOTO 23 
Z(K1)=Z(K1)+U5+U7 
IF(IP.EQ.MP) GOTO 22 

23 Z(K2)=Z(K2)+U5-U7 
GO TO 22 

20 Z(K5)=Z(K5)+U8-U9 
Z(K6)=U8+U9 
IF(IP.EQ.I) GOTO 24 
Z(K1 )=Z(K1)+U5+U7 
Z(K3)=Z(K3)+U6-U7 
IF(IP.EQ.MP) GO TO 22 

24 Z(K2)=Z(K2)+U5-U7 
Z(K4)=U6+U7 

22 CONTINUE 
JN=JN+NROW 

100 CONTINUE 

FOR MULTIPLE BORS WRITE OVER DUMMY 
SEGMENT ROWS 



C WRITE OVER ROWS AND COLS FOR FIRST DUMMY 
SEG. 

C 

DO 101 LSS=1,MT2 
LS=LSS-1 

Z(LS*NROW+NP1 1 )=(0. ,0.) 

Z(LS*NROW+NPl 1 -1 )=(0.,0.) 

101 Z(LS*NROW+NP1 1+MT)=(0.,0.) 

C 

C WRITE OVER ROWS AND COLS FOR SECOND 
DUMMY SEG. 

C 

DO 201 LSS=1,MT2 
LS=LSS-1 

Z(LS*NROW+NRF)=(0.,0.) 

Z(LS*NROW+NRF-1)s(0.,0.) 

201 Z(LS*NROW+NRF+MT)=(0.,0.) 

C 

C WRITE OVER DUMMY SEGMENTS IN THE CASE OF 
MULTIPE WIRES 
C 

IF(NWIRES.EQ.I) GO TO 292 
NSPTS=0 

DO 291 L=1 ,NWIRES-1 
NSPTS=NSPTS+NS(L) 

DO 291 1=1 ,NROW 
C 

C SET COLS TO ZERO FOR STRIPS 
C 

IN=(NSPTS-2)*NROW 

Z(IN+l)=U0 

Z(IN+kNROW)=U0 

291 CONTINUE 

292 RETURN 
END 

SUBROUTINE 

ZMATWW(N WIRES, NW1 ,NW2,XH,YH.RH,ZH.NT,XT,AT,A, 

Z) 

c 

C ••• MODS FOR DIPOLE RECTANGULAR 
COORDINATES ARE SENT OVER "* 

C COMPUTE MATRIX ELEMENTS FOR MULTIPLE 
WIRES USING GLISSON METHOD 
C NWIRES= # OF WIRES; NW1(I),NW2(I) ARE START 
AND STOP POINTS. 

C &&&&& WIRES NEED NOT BE IN THE Z=0 PLANE 
&&&&& 

C 

COMPLEX C EXP, Z(1 0000), CON, CMPLX 
COMPLEX U,U0,PSI,SUM1 ,SUM2,SUM3,U5,U6,U7 
DIMENSION 

RH(400)^H(400),XT(4),AT(4),XH(400),YH(400) 
DIMENSION D1 (400), SI (400), Cl (400) 

DIMENSION 

NW1 (4),NW2(4),NS(4),XS1 (400), YS1 (400), ZS1 (400) 
DIMENSION CU(400),SU(400) 

Pl=3. 141592 
PI2=2.*PI 
PD2-PI/2. 

U0=(0.,0.) 

U=(0.,1.) 

c 

C DEFINE GEOMETRY TERMS FOR THE WIRE. 
XH,YH,RH,ZH ARE ALL KNOWNS. 

C 
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DO 5 L=1 ,NWIRES 
5 NS(L)=NW2(L)*NW1(L)+1 
NS1=NW2(NWIRES)-NW1 (1 ) 

NPS=NS1+1 
DO 10 N=2,NPS 
N0=N-1 
|=NW1 (1)+N-1 
12 = 1*1 

XS1(N0)=.5*(XH(I)+XH(I2)) 

YS1 (N0)=.5*(YH(I)+YH(I2)) 

2S1 (N0)=.5*(ZH(I)+ZH(I2)) 

DX=XH(I)-XH(I2) 

DY=YH(I)-YH(I2) 

D1 (NO)=SORT (DX**2+DY**2) 

CU(N0)=DX/D1(N0) 

SU(N0)=DY/D1(N0) 

S 1 (N0)= 1 .0 
C1(N0)=0.0 
10 CONTINUE 
NROW=NS1-1 
JN=-1-NROW 
DO 500 JO=1,NS1 
KQ=2 

IF(J0.E0.1 ) KO=1 
IF(JO.EO.NS1 ) KO=3 
01=D1 (JO)/2. 

DO 22 IP=1,NS1 
LO=0 

IF(IP.EO JO) 10=1 

CON=Dl(IP)/4./PrU 

SUM1=U0 

SUM2=U0 

SUM3=U0 

AA=CU(IP)*CU(JO)+SU(IP)*SU(JO) 

C 

C H INTEGRATION - SUBTRACT OUT SINGULARITY IF 
JO=IP 

AND SAME SEGMENT (N1 1=N21): DESIGNATED LO=1 

DO 100 1=1, NT 
H=OrXT(l) 

IF(LO.NE.I) GO TO 40 
RP=SORT(H**2+A*A) 

GO TO 50 

40 XP=XS1(JO)+H*CU(JO)*S1(JO) 

YP=YS 1 (JO)+H*SU(JO)*S1 (JO) 

RP=SORT((XS1 (IP)-XP)**2+(YS1 (IP)-YP)**2+(ZS1 (JO)-ZS 
1(IP))—2) 

50 PSI=CEXP(CMPLX(0.,*RP)) 
c IF(LO.NE.I) GO TO 60 
c PSI=PSI*(1.,0.) 

60 PSI=PSI/RP 

SUM1=SUM1+PSI*AA*AT(l)/4. 
SUM2=SUM2+PSI*AA*AT(l)/2. Ol (JO)*H 
SUM3=SUM3+PSI*AT(I)/D1 (IP)/D1 (JO) 

100 CONTINUE 
X1=0. 

X2=0. 

X3=0. 

c IF(LO.NE.I) GOTO 200 
c DO=D1(JO)/2. 
c SDO=SORT (DO**2+A**2) 
c XIO=ALOG(SDO+DO)-ALOG(SDO-DO) 
c X1=XI0’AA/4. 



c X3=XI0/D1 (IP)/D1 (JO) 
c 200 CONTINUE 

SUM 1 =(SUMr 01 +X 1 )*CON 

SUM2=(SUM2*01+X2)*C0N 

SUM3=(SUM3*01+X3)’C0N 

U5=SUM1-SUM2 

U6=SUM1+SUM2 

U7=-SUM3 

K1=IP+JN 

K2=K1+1 

K3=K1+NROW 

K4=K2+NROW 

GOTO (18, 20, 19), KO 

18 CONTINUE 
IF(IP.EO.I) GO TO 21 
Z(K3)=Z(K3)+U6-U7 
IF(IP.EQ.NS1 ) GO TO 22 

21 Z(K4)=U6+U7 
GO TO 22 

19 CONTINUE 
IF(IP.EO.I) GOTO 23 
Z(K1 )=Z(K1)+U5+U7 
IF(IP.EO.NSI) GO TO 22 

23 Z(K2)=Z(K2)+U5-U7 
GOTO 22 

20 CONTINUE 
IF(IP.EO.I) GOTO 24 
Z(K1)=Z(K1)+U5+U7 
Z(K3)=Z(K3)+U6-U7 
IF(IP.EO.NSI) GOTO 22 

24 Z(K2)=Z(K2)+U5-U7 
Z(K4)=U6+U7 

22 CONTINUE 
JN=JN+NROW 

500 CONTINUE 

IF(NWIRES.EO.I) GO TO 297 
NSPTS=0 

DO 291 L=1 ,NWIRES*1 
NSPTS=NSPTS+NS(L) 

DO 291 1=1 ,NROW 
C 

C USE THE SIMPLIFIED METHOD OF HANDLING 
MULTIPLE WIRES 

C SET COLS TO ZERO FOR STRIPS 
C 

IN=(NSPTS-2)*NROW 

Z(IN+l)=U0 

Z(IN+l+NROW)=U0 

C 

C SET ROWS TO ZERO FOR STRIPS 
C 

IN=NSPTS-2+(l-1)*NROW 
Z(IN+1 )=U0 
Z(IN+2)=U0 
291 CONTINUE 
C 

C SET DIAGONALS TO 1 FOR STRIPS 
C 

NSPTS=0 

DO 296 L=1,NWIRES*1 
NSPTS=NSPTS+NS(L) 

I0=NSPTS*2 

M=IO*NROW+NSPTS-2 

Z(I1+1)=(1 M 0.) 

Z(l1+NROW+2)=(1 .,0.) 
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296 CONTINUE 

297 RETURN 
END 
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